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

基于遺傳算法的衛星攻擊軌道優化方法

2016-11-02 00:38:31李玉玲王藝鵬
系統工程與電子技術 2016年5期
關鍵詞:規劃

趙 琳,李玉玲,劉 源,郝 勇,王藝鵬

(哈爾濱工程大學自動化學院,黑龍江 哈爾濱 150001)

基于遺傳算法的衛星攻擊軌道優化方法

趙琳,李玉玲,劉源,郝勇,王藝鵬

(哈爾濱工程大學自動化學院,黑龍江哈爾濱 150001)

在空間攻防與衛星對抗中,當目標衛星周圍有若干小衛星以編隊形式對其繞飛時,為使攔截衛星成功擊中目標衛星,并且避開編隊小衛星的防御區,必須對攔截衛星攻擊軌道進行規劃。尋找到一條既能滿足安全性、快速性,又能節省燃料的最優路徑。而利用經典數學規劃方法,如序列二次規劃方法,雖能尋找到最優路徑,但并不適應于解決空間對抗中復雜攻防環境模型下的軌道規劃問題。為此本文提出基于遺傳算法的攔截衛星攻擊軌道尋優方法。建立目標衛星編隊小衛星的動態防御模型作為環境模型,采用可變長度實數編碼方式,根據攻擊軌道安全性、快速性、燃料消耗最少等要求建立綜合適應度函數,并對遺傳算子及置換運算方法進行設計。通過仿真驗證,本文提出的軌道優化方法能夠尋求到最優攻擊路徑,并且算法收斂速度較快。

軌道優化;遺傳算法;動態模型;攔截衛星

網址:www.sys-ele.com

0 引 言

近年來,隨著空間技術的迅猛發展,空間平臺在現代戰場中的作用越來越大,對戰局的發展甚至戰爭結果的影響越來越明顯。特別是當各種軍用衛星融入信息化作戰體系,成為這個體系的重要組成環節時,反衛星作戰已經無法避免,它成為信息化戰爭最有效、最關鍵的破門之戰。為此,各國都在迅速提升空間戰爭能力和反衛星能力。攔截衛星,作為空間攔截的主要方式之一,在攔截過程中,必須安全突破編隊小衛星的動態防御區,才能接近并攔截目標衛星。因此,其軌道規劃是成功攔截的關鍵因素。

軌道規劃[1 5]的主要任務是在具有障礙物的環境中,按照一定的評價標準,尋找一條從起始狀態到達目標狀態的無碰撞路徑。常見的路徑規劃方法有基于勢場法、構型空間法、神經網絡和可視頂點圖法等,但這些算法都存在著一些算法本身的局限性。許多算法運算量過大,實時性差[6 10]。在衛星軌道設計中,比較常用的算法通常是經典的數學規劃算法(如序列二次規劃方法),此方法比較適用于無障礙環境下的軌道規劃。其基本思想是利用二次近似在給定的近似點逐漸得到更優的迭代點,可通過求解二次規劃子問題得到,并且在當前迭代點處通過求解一系列的二次規劃子問題使迭代點逐漸靠近原始優化問題的最優解。而在空間對抗中,由于目標衛星周圍有若干伴星構成防御系統對其繞飛,因此環境模型較為復雜,在此條件下利用序列二次規劃方法很難尋找到一條能夠避開目標衛星防御區的安全路徑。而遺傳算法具有很強的全局尋優能力和適應性,被廣泛地應用于路徑規劃問題,或者用來對原有路徑規劃方法進行改進[11 13]。在利用遺傳算法解決攔截衛星軌道規劃問題中,文獻[13]是利用擴大編隊小衛星的防御半徑的方法將動態模型轉化成靜態模型,這種方法雖能滿足安全性的要求,但是這種轉化使算法搜索空間受到限制,搜索到的最優路徑可能并不能很好滿足路程最短的要求。

本文通過建立編隊小衛星在相對軌道的動學力模型,得到攔截衛星從起始時刻到達攻擊目標的時刻編隊小衛星的動態防御模型,在此動態環境模型的基礎上對遺傳算法進行設計,通過計算機仿真證明該方法能解決動態環境中的攔截衛星軌道規劃問題。

1 問題描述

目標衛星周圍有若干個小衛星以編隊形式對其繞飛,當攔截衛星攻擊目標衛星時,小衛星可通過軌道機動方式打擊攔截衛星。攔截衛星必須安全躲避編隊小衛星的防御才能成功擊中目標衛星。因此,攔截衛星應尋找到一條安全路徑,使該路徑滿足燃料消耗最少,攻擊時間最短的要求。

2 動力學模型的建立

為建立目標衛星、編隊小衛星及攔截衛星的相對軌道動力學模型,需建立如圖1所示的相對軌道參考坐標系。其中,O-XYZ為地心慣性坐標系,OX軸在赤道面內指向春分點方向,OZ軸垂直于赤道面指向北極方向,OY軸由右手法則確定。o-xyz為目標衛星軌道坐標系,原點o在目標衛星質心,x軸方向為目標衛星運動方向,y軸方向為目標衛星軌道角動量方向,z軸與x軸、y軸滿足右手定則。

圖1 相對軌道參考坐標系

在目標衛星軌道坐標系中,攔截衛星及伴星的相對動能可表示為

式中,L為拉格朗日函數;m為各星質量;r為各星相對目標衛星的位置矢量;ω為目標衛星軌道坐標系相對地心慣性坐標系轉動的角速度矢量;a0為目標衛星軌道坐標系原點的加速度矢量;V(r,ρ,t)和V(r,t)分別為伴星和目標衛星的單位質量勢能。由拉格朗日方程可知,相對動力學方程可表示為

式中,F為各星所受的不能用勢函數描述的其他攝動力和主動控制力。將式(1)帶入到式(2)中,得到在目標相對軌道坐標系o-xyz中的標量表達式為

式中,(ax,ay,az)為攔截衛星在相對軌道坐標系中的運動加速度在各軸上的分量。在二體引力作用下,當目標衛星運行在近似圓軌道,并且編隊小衛星與目標衛星之間的距離相比目標衛星軌道長半軸為小量時,編隊小衛星在相對參考坐標系中的運動可由C-W方程描述:

由C-W方程可知,若已知各小衛星在相對運動軌道中的初始狀態及單位質量作用力,就可得出其相對軌道動力學模型,進而可推算出任意時刻編隊小衛星在相對軌道中的位置和速度。初始狀態選取不同時,編隊形式也不同,當小衛星以空間圓編隊繞飛時,應滿足

當小衛星以同軌道面編隊繞飛時,應滿足

本文采用8顆衛星組成編隊小衛星,目標衛星運行軌道為圓軌道,且無軌道機動,其軌道參數如表1所示。編隊小衛星相對軌道初始狀態如表2所示。編隊小衛星在軌道中的運動模型如圖2所示。

圖2 編隊小衛星在軌道中的運動模型

表1 目標衛星軌道參數

表2 編隊小衛星相對軌道初始狀態

3 算法設計

3.1染色體可變長度編碼

本文采用可變長實數編碼方式[14],攔截衛星利用末段軌道脈沖變軌方式逐漸接近目標衛星。此編碼方式指在種群進化過程中利用交叉操作使染色體長度(即編碼長度)發生改變,

在遺傳搜索過程中逐漸趨近最優解,并且可加快算法的收斂速度[15]。在脈沖施加時刻分別為(t1,t2,…,tn),速度變化量為(Δυ1,Δυ2,…,Δυn),其中υi=(Δυix,Δυiy,Δυiz)。假設tn+1時刻攔截衛星擊中目標衛星,那么tn時刻施加脈沖時攔截衛星速度υn可由tn時刻攔截衛星位置以及目標衛星位置推算出來,故路徑編碼可表示為(t1,Δυ1,t2,Δυ2,…tn-1,Δυn-1,tn)。

在大范圍空間內全局規劃時,與傳統的二進制編碼方式、柵格序號編碼方式[16 18]以及空間直角坐標編碼方式[19]相比,本文所采用的時間與速度變化量串聯的可變長編碼方式具有編碼長度短、簡明、直觀、收斂速度快等優點,更適用于在動態障礙環境中搜尋最優路徑,并且本文可通過速度變化量直接推算出變軌所需的脈沖強度,可行性強。

3.2初始種群的產生

采用隨機生成的方式產生NUM個個體構成初始種群。每個個體記為Pj(j=1,2,…,NUM)。并且使t1<t2,…,<tn,t1=0。

3.3適應度函數的建立

3.3.1描述路徑安全性的適應度函數

為保證路徑安全,應避免路徑中某時刻攔截衛星的位置在編隊小衛星的防御范圍內,因此建立路徑安全約束適應度函數如下:

lji表示第j條路徑ti時刻到ti+1時刻的位置矢量;Φji為編隊小衛星的防御區。ηji為1時表示路徑安全,為0時表示路徑穿越防御區。由此可知,此適應度函數越大,該路徑越安全。

3.3.2描述變軌次數的適應度函數

式中,Nj表示第j條路徑攔截衛星的脈沖變軌次數。該適應度函數值越小,變軌次數越少,攔截衛星的燃料消耗越少。

3.3.3描述路徑長度的適應度函數

式中,(xji,yji,zji)表示第j條路徑ti時刻攔截衛星在相對軌道坐標系中的位置;d((xji,yji,zji),(xj(i+1),yj(i+1),zj(i+1)))表示第j條路徑ti時刻到ti+1時刻的路徑長度。該適應度函數值越小表明路徑越短。

3.3.4描述時間長度的適應度函數

式中,Tj表示第j條路徑從起始時刻到擊中目標衛星時刻所需時間。該適應度函數越小,擊中目標衛星所用時間越短。

3.3.5描述燃料消耗的適應度函數

式中,υji表示第j條路徑第i次施加脈沖時攔截衛星的速度變化量。該適應值越小,攔截衛星變軌消耗的燃料越少。

3.3.6適應度函數歸一化處理

利用適應性權重方法對上述適應度進行歸一化處理,均衡各個適應度函數值對算法收斂性的影響[20 23]。當前種群中第k個適應度函數的最大值和最小值定義如下:

則第k個適應度函數的適應性權重表示為

故得到綜合適應度函數如下:

綜上可知,該綜合適應度函數值越大,個體越優良。

3.4遺傳算子

3.4.1選擇算子

選擇算子采用輪盤賭方式,根據種群中每個個體的適應度值的比例來確定該個體被選擇的概率。這種選擇方法具有隨機采樣的特點,使種群中每個個體都有被選中的概率。首先計算個體適應度值的比例:

式中,fit(Pj)為第j個個體的適應度值。若滿足

則第j個個體被選中,其中M為[0,1]間的隨機數。

3.4.2交叉算子

交叉運算采用單點交叉方式,使2個相互配對的個體交換部分基因。具體操作為:對種群中的所有個體進行兩兩配對,分別在每個個體中隨機選擇節點作為交叉點,對每一對相互配對的個體,按照交叉概率在交叉處相互交換部分基因從而產生新的個體。圖3為2個父代經雜交操作后得到新個體的過程。

圖3 交叉操作示意圖

3.4.3變異算子

變異算子采用非均勻變異方式,具體操作為:對于給定的父代x,如果按照變異概率它的元素xk被選中進行變異,結果的后代是x=(x1,…,xk,…,xn),其中xk從下面兩種可能的方案中隨機選擇:

式中,xUk,xLk分別為個體中第k個元素所能取到的最大值和最小值;r為[0,1]中的隨機數;T為最大遺傳代數;t為當前遺傳代數;b是描述飛均勻程度的參數。該變異方式可使算子在早期均勻第搜索解空間,而到了晚期則在很小的區域進行搜索,可加快算法的收斂速度。

3.4.4置換運算

在遺傳算法中種群通過進化可得到越來越多的優良個體,但是由于遺傳算子操作中的隨機性,可能導致某些優良個體被破壞,因此可通過置換操作,將種群中的較優個體保留,加快算法的收斂速度。具體操作為:首先將父代個體x(t)進行選擇、交叉、變異操作后得到子代個體x(t),根據置換比例Pm從父代種群中選擇M×Pm個最優個體,再從子代種群中選擇M×(1-Pm)個個體共同組成新種群x(t+1)。

3.5遺傳算法過程

算法流程如圖4所示。

圖4 算法流程圖

步驟1種群初始化。

步驟2計算個體適應值。

步驟3判斷是否滿足算法終止條件:如果不滿足,程序繼續;如果滿足,程序結束。

步驟4將父代個體進行選擇、交叉、變異操作,得到子代個體。

步驟5根據置換比例從父代種群中選擇M×Pm個最優個體,再從子代種群中選擇M×(1-Pm)個個體共同組成新種群,返回步驟2。

4 仿真校驗

本文以Matlab為實驗平臺,以8顆小衛星組成的防御模型為例,攻擊衛星在相對軌道中的起點(單位:m)設為(6 000,10 000,10 000),目標衛星位置為(0,0,0),種群規模為200,交叉概率為0.5,變異概率為0.02,置換比例為0.02,迭代次數為1 000代。為驗證遺傳算法在復雜環境的攻防對抗中對軌道規劃的有效性及優越性,本文將遺傳算法的搜索結果與序列二次規劃方法相比較。利用遺傳算法的仿真結果如圖5~圖7所示。利用序列二次規劃方法仿真結果如圖8和圖9所示。

圖5 遺傳算法迭代1 000次的最優路徑

圖6 遺傳算法得到的最優路徑中各小衛星與攔截衛星相對距離減去防御半徑后隨時間的變化

圖7 隨遺傳算法迭代次數增加適應值變化過程

經仿真驗證,第425 s時攔截衛星擊中目標衛星。表3為經過Matlab仿真得到的遺傳算法軌道規劃的脈沖施加時間及大小。圖5為遺傳算法迭代1 000次后得到的最優路徑,圖中陰影部分攔截衛星經過各小衛星附近時的危險防御區域。從圖5中可以看出,該路徑中沒有穿越編隊小衛星防御區的點,滿足路徑安全性要求。路徑曲線較為平滑,每次施加脈沖時,攔截衛星變軌角度較小,與以路徑點的坐標為編碼方式尋求到的最優路徑相比,更能到達滿足節省燃料的要求。圖6為最優路徑中各編隊小衛星與攔截衛星的相對距離減去小衛星防御半徑后隨時間變化的過程,由圖6可知,攔截衛星安全避開了編隊小衛星的防御區,沒有發生碰撞,并且由圖示縱坐標可看出攔截衛星從起點到擊中目標所用時間為419 s,可滿足攻擊時間最短的要求。圖7為進化過程中種群適應值變化,算法進化到第150代左右尋找到最優路徑,可見算法收斂性和穩定性較好。

表3 遺傳算法軌道規劃的脈沖施加時間及大小

與文獻[13]結果相比,本文搜索到的最優路徑變軌次數少,且變軌角度小,更能到達節省燃料的要求。并且與文獻[13]中的編碼方式相比較,本文采用的編碼方式,可直接得出攔截衛星攻擊所用時長及每次變軌的時刻,可行性更強。本文采用的編隊小衛星的動態防御模型使算法擁有更大的搜索空間,搜索到的最優路徑更可靠。

圖8 序列二次規劃方法得到的最優路徑

圖9 序列二次方法得到最優路徑中各小衛星與攔截衛星相對距離減去防御半徑后隨時間的變化

圖8中陰影部分攔截衛星經過各小衛星附近時的危險防御區域。由圖8和圖9可知,當目標衛星周圍有若干伴星組成防御系統對其繞飛時,利用序列二次規劃方法尋求到攻擊路徑有穿越編隊小衛星防御區的現象,雖然能做到能量最優,但并不能滿足安全性方面的要求,因此在復雜環境下的攻防對抗中,序列二次規劃方法仍然存在一定局限性。

5 狀態誤差分析

當攔截衛星在相對軌道坐標系中未施加脈沖作用時,式(4)可簡化為

利用式(4)及式(18)就可得到利用脈沖加速度進行軌道機動與未施加脈沖時的誤差模型:

利用式(19)以及表3中的仿真數據,就可分析到達每個目標點的狀態誤差,如表4所示。

表4 每個目標點的狀態誤差

6 結 論

針對空間攻防對抗中的攻擊衛星軌道規劃問題,本文基于遺傳算法,在編隊小衛星的動態防御模型的基礎上,提出了攔截衛星尋求最優攻擊路徑的方法。采用可變長實數編碼方式對路徑進行編碼,根據安全性、燃料最省和攻擊時間最短等約束條件建立適應度函數,并且對遺傳算子和置換運算方法進行了設計。仿真結果表明,本文提出的算法能夠尋求到既滿足安全性,同時滿足燃料最省和攻擊時間最短等要求的最優路徑,與其他路徑規劃方法相比,本文提出的算法更適合在復雜環境模型下進行軌道規劃。

[1]Meng X Q,Zhao Y N,Xue Q,Application of genetic algorithm in path planning[J].Computer Engineering,2008,34(16):215 220.(孟憲權,趙英男,薛青.遺傳算法在路徑規劃中的應用[J].計算機工程,2008,34(16):215-220.)

[2]Bulut O,Tasgetiren M F.An artificial bee colony algorithm for the economic lot scheduling problem[J].International Journal of Production Research,2014,52(4):74-76.

[3]Al-Dallal A,Rasha S A,El-Haddadeh R.IR with and without GA:study the effectiveness of the developed fitness function on the two suggested approaches[J].International Journal of Applied Metaheuristic Computing,2013,4(1):1-20.

[4]Liu T F,Cheng R Y.Mobile robot path planning based on genetic algorithm[J].Computer Engineering,2008,34(17):214 215.(劉天孚,程如意.基于遺傳算法的移動機器人路徑規劃[J].計算機工程,2008,34(17):214-215.)

[5]Zheng Q,Sha J X,Shu H,et al.A variant constrained genetic algorithm for solving conditional nonlinear optimal perturbations[J].Adυances in Atmospheric Sciences,2014,3(1):219-229.

[6]Wang C,Gao J H.A differential evolution algorithm with cooperative coevolutionary selection operation for high-dimensional optimization[J].Optimization Letters,2014,82(2):477-492.

[7]Keshavarz S,Khoei A R,Molaeinia Z.Genetic algorithm-based numerical optimization of powder compaction process with temperature-dependent cap plasticity model[J].The International Journal of Adυanced Manufacturing Technology,2013,64(5):1057-1072.

[8]Ren Y,Cui P Y,Luan E J.Low-thrust trajectory optimization based on annealing-genetic algorithm[J].Journal of Astronautics,2007,28(1):162-166.(任遠,崔平遠,欒恩杰.基于退火遺傳算法的小推力軌道優化問題研究[J].宇航學報,2007,28(1):162-166.)

[9]Wang R L,Okazaki K.An improved genetic algorithm with conditional genetic operators and its application to set-covering problem[J].Soft Computing,2013,11(7):687-694.

[10]Lu T,Zhu J.A genetic algorithm for finding a path subject to two constraints[J].Applied Soft Computing Journal,2013,13(2):115-121.

[11]Zhang Y,Gong D W,Sun X Y,et al.Adaptive bare-bones particle swarm optimization algorithm and its convergence analysis[J]. Soft Computing,2014,18(7):1337-1352.

[12]Sitansh S,Harjinder S.Genetic algorithm optimization of laser pulses for molecular quantum state excitation.[J].The Journal of Chemical Physics,2010,132(6):64-108.

[13]Deng H,Zhong W C,Sun Z W,et al.Method resrarch of satellite attacking path planning based on genetic algorithm[J].Journal of Astronautics,2008,30(4):1587-1592.(鄧泓,仲惟超,孫兆偉,等.基于遺傳算法的衛星攻擊路徑規劃方法研究[J].宇航學報,2008,30(4):1587-1592.)

[14]Taherdangkoo M,Paziresh M,Yazdi M.An efficient algorithm for function optimization:modified stem cells algorithm[J].Central European Journal of Engineering,2013,3(1):36-50.

[15]Song X M,Zhan C S,Xia J.Integration of a statistical emulator approach with the SCE-UA method for parameter optimization of a hydrological model[J].Chinese Science Bulletin,2012,57(26):3397-3403.

[16]Leung C S,Lam P M,Situ WC.A graphics processing unit accelerated genetic algorithm for affine invariant matching ofbroken contours[J].Journal of Signal Processing Systems,2012,66(2):105-111.

[17]Cao X B,He D L.Error analysis of Hill equation-based formation initialization[J].Flight Dynamics,2008,26(6):84-88.(曹喜濱,賀東雷.基于Hill方程的編隊初始化誤差分析[J].飛行力學,2008,26(6):84-88.)

[18]You S H,Lee Y H,Lee W J.Parameter estimations of a storm surge model using a genetic algorithm[J].Natural Hazards,2012,60(3):1157-1165.

[19]Shan H B,Zhou S H,Sun Z H.Research on assembly sequence planning based on genetic simulated annealing algorithm and ant colony optimization algorithm[J].Assembly Automation,2009,29(3):249-256.

[20]K?ker R.A neuro-simulated annealing approach to the inverse kinematics solution of redundant robotic manipulators[J].Engineering with Computers,2013,29(4):507-515.

[21]Zhang Q,Xie Z P,Liu Z T.Minmal reduct computing based on variable length coding genetic algorithm[J].Mini-Micro System,2001,22(9):1055-1057.(張卿,謝志鵬,劉宗田.基于變長編碼遺傳算法的最小縮減計算[J].小型微型計算機系統,2001,22(9):1055-1057.)

[22]Liu Y,Shen Y,Xing L,et al.Equal life design method of operationally responsive on-board electronic systems[J].Acta Aeronautica et Astronautica Sinica,2014,35(6):1673-1683.(劉源,沈毅,邢雷,等.快速響應衛星電子系統等壽命設計方法[J].航空學報,2014,35(6):1673-1683.)

[23]Lozano P T,Wesley M A.An algorithm for planning collisionfree paths among polyhedral obstacles[J].Communications of the ACM,1979,22(10):560-570.

劉源(1984-),通訊作者,男,講師,博士,主要研究方向為空間攻防、星載電子系統、飛行器任務規劃。

E-mail:spacead@163.com

郝勇(1979-),男,講師,博士,主要研究方向為集群航天器的姿態協同控制。

E-mail:haoyong@hrbeu.edu.cn

王藝鵬(1990-),男,博士,主要研究方向為神經網絡與組合導航系統。

E-mail:wangyipeng@hrbeu.edu.cn

Optimization method research of satellite attaching track planning based on genetic algorithm

ZHAO Lin,LI Yu-ling,LIU Yuan,HAO Yong,WANG Yi-peng
(College of Automation,Harbin Engineering Uniυersity,Harbin 150001,China)

In the warfare of anti-satellite,when formation micro-satellite is flying around the target satellite,in order to avoid the defensive area of micro-satellite,it is needed to find a track security track which needs less fuel.A method to find the best attacking track based on the genetic algorithm is proposed.The dynamic defensive model of the target micro-satellite is built as the environment model,using variable length encoding real numbers.Establish the comprehensive fitness function based on security,rapidity and fuel,and design the selection operator,crossover operator,mutation operator and replacement calculation method.Numerical simulations show that the proposed method can get the best attacking track,and have a faster convergence rate.

track optimization;genetic algorithm;dynamic model;attacking satellite

V 19

A

10.3969/j.issn.1001-506X.2016.05.22

1001-506X(2016)05-1114-07

2015-04-30;

2015-07-24;網絡優先出版日期:2015-12-29。

網絡優先出版地址:http://www.cnki.net/kcms/detail/11.2422.TN.20151229.1748.010.html

黑龍江省博士后科研啟動金(LBH-Q14054)資助課題

趙琳(1968-),男,教授,博士,主要研究方向為衛星導航、組合導航。

E-mail:zhaolin@hrbeu.edu.cn

李玉玲(1992-),女,碩士研究生,主要研究方向為控制科學與工程。

E-mail:liyuling@hrbeu.edu.cn

猜你喜歡
規劃
我們的規劃與設計,正從新出發!
房地產導刊(2021年6期)2021-07-22 09:12:46
“十四五”規劃開門紅
“十四五”規劃建議解讀
發揮人大在五年規劃編制中的積極作用
規劃計劃
規劃引領把握未來
快遞業十三五規劃發布
商周刊(2017年5期)2017-08-22 03:35:26
基于蟻群算法的3D打印批次規劃
多管齊下落實規劃
中國衛生(2016年2期)2016-11-12 13:22:16
十三五規劃
華東科技(2016年10期)2016-11-11 06:17:41
主站蜘蛛池模板: 国产日本欧美在线观看| 亚洲人成网址| 67194成是人免费无码| 久久特级毛片| 色哟哟国产精品| 亚洲精品中文字幕无乱码| 无码中文AⅤ在线观看| 欧美第一页在线| 久久伊人久久亚洲综合| 欧美中日韩在线| 国产精品久久久久久影院| 青青青国产免费线在| 国产二级毛片| 在线观看免费国产| 国产素人在线| 波多野结衣一区二区三区四区视频 | 国产天天色| 亚洲欧洲日韩综合色天使| 麻豆a级片| 久久熟女AV| 亚洲最大综合网| 91免费国产高清观看| 国产精品福利在线观看无码卡| 日韩精品久久无码中文字幕色欲| 亚洲一区二区成人| 日韩天堂网| 久久香蕉欧美精品| 狠狠干欧美| 伊人网址在线| 国产精品久久久精品三级| 波多野结衣无码中文字幕在线观看一区二区| 激情爆乳一区二区| 久久亚洲欧美综合| 国产欧美日韩精品综合在线| 青青青草国产| 精品国产成人三级在线观看| 午夜高清国产拍精品| 国产福利一区在线| 久久精品亚洲热综合一区二区| 园内精品自拍视频在线播放| 国产成年无码AⅤ片在线| 欧美国产精品不卡在线观看| 男女猛烈无遮挡午夜视频| 亚洲国产成人久久77| 免费人成又黄又爽的视频网站| 一本大道无码高清| 久久亚洲国产视频| 毛片三级在线观看| 免费在线国产一区二区三区精品| 免费全部高H视频无码无遮掩| 国产成人调教在线视频| 久久中文电影| 99国产精品国产| 红杏AV在线无码| 国产99热| 91在线精品免费免费播放| 老司机久久精品视频| 有专无码视频| 欧美日本在线一区二区三区| 亚洲无码高清免费视频亚洲 | 国产乱子伦一区二区=| 欧美成人看片一区二区三区 | 114级毛片免费观看| 国产精品视屏| 亚洲视频四区| 欧美精品亚洲精品日韩专区| 亚洲AⅤ永久无码精品毛片| 高清精品美女在线播放| 色婷婷狠狠干| 亚洲浓毛av| 亚洲中文精品人人永久免费| 欧美午夜理伦三级在线观看| 无码啪啪精品天堂浪潮av| 国产日韩欧美成人| 日韩专区欧美| 精品人妻无码中字系列| 18禁色诱爆乳网站| 国产欧美成人不卡视频| 动漫精品啪啪一区二区三区| 国产免费一级精品视频| 国产麻豆精品久久一二三| 色窝窝免费一区二区三区|