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

基于多目標帝國競爭算法的進場排序與調度

2021-07-07 11:36:52張軍峰游錄寶楊春葦胡榮
航空學報 2021年6期
關鍵詞:排序優化

張軍峰,游錄寶,楊春葦,胡榮

南京航空航天大學 民航學院,南京 210016

空中交通需求的持續增長與可用空域資源的長期受限,給空中交通管理(Air Traffic Management,ATM)帶來了新的機遇與挑戰。以中國繁忙機場的終端區運行為例,完全依賴于管制員指令引導的方式,容易導致管制工作負荷激增、機場運行效率降低、環境影響問題突出。因此,如何有效地優化和調度時空資源,成為空中交通管理領域的研究熱點,而進場排序與調度是該領域的典型問題。

進場排序與調度旨在不違反安全間隔的條件下,結合運行約束,合理高效地為進場航空器分配著陸跑道,提供最優著陸次序與時間,以期達到提升跑道容量、減少延誤、緩解管制工作負荷的目的[1]。歷經近30年,進場排序與調度的研究涵蓋了靜態優化[2]與動態優化[3],確定型優化[4]與隨機型優化[5]、單階段優化[6](直接優化著陸時間)與兩階段優化[7-8](先確定著陸次序再優化著陸時間)、是否提供相應管制建議[9](即優化時間的可達性)等方面,其求解工具或算法包括:CPLEX[1,7]、分支定界[5]、動態規劃[6]、模擬退火[10]、遺傳算法[11]、粒子群算法[12]等。

對于進場排序與調度問題,各利益相關方秉持不同的訴求:空管立足運行安全,航司著眼效率優先,機場注重容量增強,民眾關切準點運行與環境影響。因此,近年來進場排序與調度的研究重點由單目標優化逐步轉向多目標優化。Sam等[13]研究了具有不同目標函數的進場排序與調度問題,并考察了不同性能指標的解決方案之間的差異。Hong等[14]以總飛行時間和序列變化次數最小化為目標,基于混合整數線性規劃對著陸次序和時間進行優化。Zhang等[15]采用基于帕累托支配的多目標模擬退火算法(MOSA),解決了在考慮連續航班運行的條件下,最大化跑道運行能力和最小化飛行延誤的進場排序問題。帶精英策略的非支配排序遺傳算法(NSGA-II)在多目標進場排序與調度研究中應用廣泛:Mokhtarimousavi等[16]考慮了跑道容量、地面成本和環境成本,基于NSGA-II實現進場排序與調度的多目標優化;馬園園等[17]面向延誤時間、延誤架次和運行公平性的目標,利用NSGA-II求解多目標進場排序與優化;劉繼新等[18]在綜合考慮空管、航空公司和機場三方的利益需求下,提出一種時隙交換方法,并基于NSGA-II解決不同交通密度情況下的進場排序與調度問題。

上述多目標優化兼顧了諸多利益相關方的訴求,也豐富了進場排序與調度研究的內涵。然而,仍存在兩方面問題亟待解決:① 諸多優化目標之間是否存在冗余,優化目標選擇能否找到依據;② 如 何評價多目標優化的解集,能否在常用的NSGA-II之外尋求更為精準和高效的多目標優化求解算法。本文擬通過將進場排序與調度問題等效于機器調度問題,借鑒機器調度領域的研究成果,為多目標進場排序與調度的目標選擇提供依據;進一步,引入非支配排序,基于近年來廣泛應用的帝國競爭算法[19],設計多目標帝國競爭算法(Multi-Objective Imperialist Competitive Algorithm,MOICA),解決進場排序與調度問題;最后,基于通用數據集與實際運行數據構建仿真驗證的場景,一方面考慮帕累托最優解集的評價指標,對比分析本文提出的MOICA與進場排序與調度領域廣泛使用的NSGA-II及MOSA,另一方面考察了MOICA在實際管制運行中的優化效能。

1 多目標進場調度模型構建

1.1 問題等效

進場航空器的排序與調度問題可視為機器調度問題[20-21]:跑道為機器,進場航空器為工件,其進入終端空域的時間等效于工件的出現時間,航空器的最早、最晚以及計劃到達時刻,分別與工件提交時間、最后期限以及交付時間一一對應,而尾流間隔約束則與順序決定的準備時間相當。為后續闡述方便,相關變量如表1所示。

表1 變量定義

1.2 目標選擇

最小化最大完成時間min(maxCi),表明最后優化一架航空器盡早著陸,意味著跑道運行容量的最大化,該目標是機場的關注點。

最小化最大流動時間min(maxFi),表明最小化進場航空器在終端空域飛行的最長時間,意味著對各個航空公司的航空器實施公平的調度。

最小化最大延誤時間min(maxTi),表明最小化進場航空器的最大延誤,體現對航空器實施公平調度。

Cmax(S*{max})=max(C(S*{Σ}))

證明如下:

首先,設Cmax(S*{max})>max(C(S*{Σ}))。因為Cmax(S*{max})=max(C(S)),定義S為任意一個調度,而S*{Σ}為S的其中一種情況,所以假設Cmax(S*{max})>max(C(S*{Σ}))不成立。

最后,證得推論成立。

1.3 模型構建

基于上述分析,多目標進場優化排序與調度問題可以描述為:對于給定n架進場航空器,當滿足著陸時間窗約束、尾流間隔約束、著陸次序的最大位置轉換約束等,以最小化進場航空器終端區內的總飛行時間、最小化最大飛行時間、最小化進場航空器的總延誤時間為目標,優化進場航空器的著陸序列及時間。構建數學模型如下:

(1a)

(1b)

min maxFi

(1c)

(1d)

(1e)

(1f)

(1g)

(1h)

(1i)

kCPS≤K

(1j)

其中:式(1d)為進場航空器著陸時間的時間窗限制;式(1e)表示著陸次序;式(1f)確保了進場航空器前后機之間的尾流間隔約束;式(1g)與式(1h)表示進場航空器的延誤時間定義及約束;式(1i)表示進場航空器在終端空域的飛行時間;式(1j)表示著陸次序調配需要滿足的最大位置轉換約束。

2 多目標進場調度算法設計

2.1 帝國競爭算法

帝國競爭算法(Imperialist Competitive Algorithm, ICA)[19]是一種社會政治進化算法,該算法將初始解集視為若干個帝國(Emp.),每個帝國由一個殖民國家(Imp.)與若干個殖民地(Col.)組成,其中殖民國家代表該帝國內最優解,通過成本值(Costi=f(Countryi))確定。對于最小化問題,其算法簡要流程如下所示:

步驟1初始化。隨機產生初始國家集合N,根據成本值選取殖民國家,分配殖民地,組建帝國。

步驟2同化。帝國中的殖民地在一定概率下向殖民國家靠攏。

步驟3革命。帝國中的殖民國家及殖民地在一定概率下發生突變。

步驟4帝國內競爭。更新帝國內殖民國家與殖民地的成本值,成本值最小的成為新的殖民國家。

步驟5帝國間競爭。更新各帝國勢力值,將最弱帝國的最弱殖民地轉移到最強帝國中。

步驟6帝國消亡。若某帝國中所有殖民地都被掠奪,則此殖民國家淪為殖民地,該帝國消亡。

步驟7驗證終止條件。若滿足,輸出最優殖民國家,否則重復步驟2~步驟6,直到滿足終止條件。

帝國競爭算法的優點在于:通過帝國內競爭和帝國間競爭,加強深度搜索和廣度搜索,從而提升了鄰域搜索和全局優化的能力。

2.2 多目標帝國競爭算法

將帝國競爭算法引入進場排序與調度的多目標優化,相較于傳統帝國競爭算法,其特點在于:① 解集的初始化與更新策略;② 面向多目標優化的成本函數的構建;③ 面向多目標優化的帝國勢力的衡量。假設多目標帝國競爭算法(MOICA)的最大迭代次數imax,種群數量為npop,殖民國家數量nimp,MOICA的算法流程圖如圖1所示。

圖1 多目標帝國競爭算法流程示意圖

1) 初始化帝國

對于進場航空器排序與調度問題,各航空器的優化著陸時間構成一個初始解,即一個國家。生成初始解集時,需要滿足時間窗約束(ri,di),本文擬采用線性插值實現解集初始化,如式(2)所示。

Ci=ξdi+(1-ξ)rii∈I

(2)

式中:ξ∈[0,1]。

由于面向多目標優化,常規的基于目標函數構建成本函數的方式并不可行。鑒于此,本文采用非支配排序與擁擠距離實現成本函數的構建,具體包括:

① 基于初始化的解集計算目標函數值,并針對所得的目標值進行快速非支配排序。

(3)

③ 初始化殖民國家,優先從帕累托前沿(Pareto Front,PF)中選取殖民國家,若PF中解的個數小于設定的帝國個數,則從下一層級中依次選擇成本值小的解作為補充。

④ 初始化殖民地,采用輪盤賭的方式,利用式(4)將殖民地分配給殖民國家。

(4)

式中:α為選擇系數。

如此完成了初始化帝國的工作,每個帝國,即進場排序與調度的初始解集的子集,包含一個PF解和若干個初始解。

2) 同 化

同化過程是帝國內殖民地向殖民國家靠攏的過程,也即進場排序與調度的可行解向帕累托最優解學習的過程。本文采用差分進化的方式,如式(5)所示,實現帝國內的同化。

Col.←Col.+β·rand·(Imp.-Col.)

(5)

式中:β表示差分進化系數;rand為隨機數。

3) 革 命

為防止迭代過程中陷入局部最優解,對殖民地進行革命操作,文中主要采用單點變異、兩點交換和子段逆序3種算子進行革命。

單點變異算子中,隨機選擇一個進場航空器,按式(2)分配新的優化著陸時間;兩點交換算子中,隨機交換兩個航空器的優化著陸時間;子段逆序算子中,隨機選擇一部分航空器的優化降落序列,并逆序。值得注意的是,在實施上述操作需要考慮位置轉換約束(Constrained Position Shift, CPS)以及相應的尾流間隔,以確保革命后的解仍是可行解。

4) 帝國內競爭

基于同化和革命的解集更新,計算目標函數值,再次進行快速非支配排序操作,利用式(3)更新國家的成本值。進一步,比較帝國內殖民國家與殖民地的成本,若殖民地的成本優于對應的殖民國家時,那么優秀殖民地成為新的殖民國家,原來殖民國家淪為殖民地。

5) 帝國間競爭

帝國間競爭是各個帝國間殖民地再分配的過程,勢力弱小的帝國將被勢力強大的帝國逐步吞并,直至消亡。帝國勢力值計算過程如下所述:

(6)

Empj.Cost

(7)

(8)

在此過程中,式(6)綜合帝國內殖民國家與殖民地的成本計算帝國總成本,其中:μ表示勢力值系數; |Empj|表示帝國中殖民地的數目;隨后,按照式(7)和式(8)實現勢力值的量綱統一,其中式(7)的作用在于防止勢力為零的情況,本文λ取值為1.2。

2.3 多目標優化解的評價

如何評價多目標進場排序與調度優化解集的優劣,選擇兼顧各利益相關方訴求的優化方案,也是值得深入研究的問題。本文擬從帕累托最優解的收斂性和多樣性兩個角度出發,通過如下4個指標[23]評估本文提出的多目標帝國競爭算法。

1) 覆蓋率 (C-Metric, CM)

覆蓋率是衡量帕累托解集收斂性的指標,即

CCM(PFA,PFB)=

(9)

式中:PFA、PFB分別為兩組Pareto最優解集。分子表示解集B受解集A支配的解的個數,分母表示解集B中所有解個數。若CCM(PFA,PFB)=1,說明解集B完全受解集A所支配,即解集B收斂性比A差。

2) 空間分布 (SPacing, SP)

空間分布是評估帕累托解集廣泛程度的指標,即

(10)

3) 超體積 (Hyper Volume, HV)

超體積是評估帕累托解集收斂性和多樣性的綜合指標,即

(11)

式中:δ(·)表示勒貝格測度,用來測量體積;Vk表示第k個點與參考點圍成的區域超體積。超體積值越大,表明解集越好。

4) 平均理想距離(Mean Ideal Distance, MID)

平均理想距離衡量帕累托解集和理想點之間距離的指標,即

CMID(PF)=

(12)

3 仿真驗證

仿真驗證從2個角度展開,一方面利用進場排序與調度研究領域常用的基準案例——OR數據集(http:∥people.brunel.ac.uk/~mastjjb/jeb/orlib/ airlandinfo.html),考察MOICA的效果,并與多目標進場排序與調度常用的NSGA-II及MOSA進行對比;另一方面采用實際運行數據,驗證MOICA的適用性。

MOICA通過MATLAB2014b編程實現,運行環境為:Windows10操作系統,Corei7-7700處理器、3.6 GHz主頻和8 GB內存的微機。

3.1 OR數據驗證

3.1.1 運行場景

OR數據集中關于進場排序與調度共有13組案例,航空器架次從10~500架不等。由于第1~8組案例不滿足進場排序與調度的尾流間隔三角不等式,故采用第9~13組數據,分別有100/150/200/250/500架次航空器,第9~12組數據,構成大規模數據集;拆分第13組數據,構成小規模數據集,拆分的依據在于——進場航空器進入終端空域的時間具有波次效應。進一步,為了反映進場航空器的交通需求特性,引入壓縮因子Kd,即

(13)

具體的數據集信息如表2所示,數據集中airland #13.1參數信息詳情如圖2所示,其中藍色方框表示航空器進入終端空域時間,黑色星號表示著陸時間窗,紅色圓圈表示預計著陸時間。

圖2 仿真場景示意圖

表2 仿真數據集信息

實驗中:MOICA應用于大規模數據集時,imax=250、npop=100、nimp=7;應用小規模數據集時,imax=150、npop=75、nimp=5;革命概率選為0.35,選擇系數為0.9,同化系數為2,勢力值系數μ=0.2。NSGA-II中交叉概率為0.7,變異概率為0.02。MOSA中初始溫度值為1 000,冷卻系數0.98。

3.1.2 結果分析

首先,針對所有數據集進行MOICA、NSGA-II和MOSA的仿真驗證,同時利用多目標優化解的評價指標進行對比,具體如表3所示。由于元啟發式算法(如模擬退火、遺傳算法、帝國競爭算法)的優化結果與參數以及初始化息息相關,表3的結果基于20次獨立仿真的均值進行對比。

由表3可知:對于14組仿真數據集,覆蓋率指標(CCM),MOICA全面優于NSGA-II與MOSA;空間分布指標(CSP),MOICA最優,MOSA較劣;超體積指標(CHV),MOICA在12組數據集中占優;平均理想距離(CMID),MOICA在13組數據集中占優。通過計算各組數據在各類指標下的均值可見,MOICA的帕累托最優解相對于NSGA-II與MOSA而言,解集更占支配地位(CCM)、解集的分布更均勻(CSP)、解集的收斂性更好(CHV)、解集的質量更高(CMID)。進一步分析MOICA解集非占優的情況發生在小規模數據集OR#13.1/OR#13.4/OR#13.7,原因在于上述數據集的壓縮因子Kd較大,對應進場需求低,排序與調度難度較小。此時,NSGA-II與MOSA也能獲得較優解。對于最復雜數據集OR#13.6,本文提出的MOICA表現更好。

表3 OR數據集仿真結果對比

其次,選取數據組OR #9,應用MOICA、NSGA-II與MOSA分別進行500次獨立實驗。針對仿真結果計算指標:空間分布、超體積、平均理想距離,并將500組評價指標以四分位圖的方式顯示仿真結果,如圖3所示。

圖3 多目標優化解評價指標的四分位圖

由圖3可知:OR #9數據的500次獨立實驗結果進一步驗證了相對于NSGA-II與MOSA而言,MOICA的優越性。覆蓋率指標未列入,原因在于MOICA帕累托解集均占支配地位。

最后,為了考察多目標優化算法的運算效率,圖4提供了不同數據集的CPU時間示意圖,無論是本文提出的MOICA,還是對比的NSGA-II,小規模數據集的迭代次數設為150次,大規模數據集的迭代次數設為250次。

圖4 多目標優化運行時間對比

由圖4可知:一方面,本文提出的MOICA其算法效率均優于NSGA-II與MOSA;另一方面,多目標優化算法的運行效率主要依賴于迭代次數和航空器架次,當迭代次數一致時(如大規模數據集,imax=250),數據集規模越大其優化時間越長。

3.2 實際數據優化分析

3.2.1 仿真場景構建

以長沙黃花機場(四字碼ZGHA))進場運行為例,圖5顯示了2019年12月26日8:00—2019年12月27日8:00共計233條進場航空器的綜合航跡,長沙機場共有5個進港點(BEMTA, LLC, LIG, OVTAN, DAPRO)。

圖5 長沙黃花機場進場運行航跡示意圖

圖6 進場航空器飛行時間四分位圖

表4 尾流間隔標準[22]

3.2.2 仿真結果分析

利用本文提出的MOICA實現長沙機場的進場排序與調度優化,選取23:00—00:30內的21架進場航空器作為優化對象,相關參數的設定同3.1節的大規模數據集。

在本次的實例仿真中,尾流間隔分別選取表4 對應標準的1.2倍、1.5倍和1.8倍實施驗證,其總延誤時間,總飛行時間和最大飛行時間的優化結果如表5所示。

由表5可知,利用MOICA實現進場的排序與調度可以有效地提升運行效率與效益,兼顧空管、航司、機場、民眾的不同訴求。即便是將尾流間隔限制擴大到現行標準的1.8倍,相對于管制實際運行結果,總延誤時間可以降低41.2%、總飛行時間可以減少11.4%、最大飛行時間可以縮短8.6%。

表5 實際案例仿真驗證結果

進一步,為考察實際著陸序列與優化著陸序列的變化情況,圖7以時序圖的形式,展示了1.5倍標準間隔設置下,各進場航空器在不同進港點(Entry-FIX)的出現時間(黑色圓圈所示)與序列、各進場航空器的預計著陸時間(Estimate Time of Arrival, ETA)(綠色三角所示)與序列、各進場航空器的優化著陸時間(Schedule Time of Arrival, STA)(藍色菱形所示)與序列、各進場航空器的實際著陸時間(Actual Time of Arrival, ATA)(紅色方框所示)與序列。

由圖7可知,1.5倍的間隔標準與實際管制運行時的安全間隔相仿,引入MOICA的優勢在于通過調整著陸次序,最大可能地逼近設置的尾流安全間隔,從而實現總延誤時間、總飛行時間以及最大飛行時間地減少。至于管制員在實際運行中未能應用最佳著陸序列的原因,可以通過不同時間點進場航空器態勢加以解釋,如圖8所示。

圖8展示了從23:30—23:55這段時間內,10架進場航空器的運行態勢。航空器的選擇是基于圖7的分析。

圖8 進場航空器不同時刻的態勢示意圖

由圖7可知,進場航空器可以粗略分為3個“波次”,第1波共有12架;第2波共有4架;第3波也是4架。前12架中,可以忽略第一架與最后一架著陸的航空器(均由DAPRO進港),因此余下10架的構成為:2架次DAPRO進港,1架次BEMTA進港,2架次LLC進港,3架次LIG進港,以及2架次OVT進港。管制員對于西向進港(LLC/BEMTA)航空器的調配策略是等待;對于北向(DAPRO)和南向(LIG)進港航空器使用向東側偏離策略;至于東向進港(OVTAN)航空器則向北側調配。23:35和23:40時,管制員發現著陸需求較為旺盛,于是分別按既定策略引導進場航空器作機動飛行。此時,管制員引導主要依賴經驗,兼之飛行員具體操控航空器時對管制員指令的依從度無法預知,最終導致無法充分利用終端空域以及跑道的容量。因此,MOICA可以為管制員提供著陸序列與時間的決策支持,至于如何給定確切的管制指令,本文擬借助于點融合系統 (Point Merge System, PMS) 實現。

圖7 進場航空器著陸時序對比示意圖

基于統計分析獲取飛行時間或者采用四維航跡預測[24],能夠實現進場運行“可見”目標;采用本文提出的MOICA兼顧各利益相關方的訴求開展排序與調度,可以實現進場運行“可優”目標;如何基于優化的著陸序列與時間向管制員提供決策支持,方能實現進場運行“可達”目標。為此,本文擬通過PMS解決優化結果“可達”的問題。圖9提供了引入PMS后的進場程序,以及基于該程序和1.2倍間隔標準優化結果的進場航空器運行示意圖(注:1.5倍間隔標準的配備需要通過引入標準等待實現)。

由圖9可見,PMS的引入具備如下優勢:①進 場排序與調度優化結果的可達,根據優化著陸序列與時間反推PMS程序結構下的航空器四維軌跡,為管制員提供決策建議;② 管制決策建議簡單有效,將傳統的管制指令簡化為排序支路上的轉彎,降低了管制工作負荷;③ 由于管制引導僅發生在排序支路上,提升了飛行員的態勢感知與情境意識;④ 飛行軌跡更為規整有序,有效降低了終端空域進場運行的復雜性。值得注意的是,由于引入PMS導致預計到達時間與著陸時間窗變化,因此表5的相關結論需要另行計算。

圖9 引入PMS的進場排序與調度運行示意圖

4 結 論

1) 借鑒機器調度領域的研究成果,梳理和刪減進場排序與調度問題所對應的目標函數,最終確定3個指標:最小化總延誤時間、最小化總飛行時間、最小化最大飛行時間。從理論角度分析和證明,上述目標能夠同時滿足空管、航司、機場以及民眾的不同訴求。

2) 面向3個目標,考慮著陸時間窗約束、尾流間隔約束、航空器著陸次序的最大位置轉換約束等,構建進場排序與調度多目標優化模型,以實現管制工作負荷降低、跑道運行容量增強、航空器正點率提升、各航司調度公平的目的。

3) 基于非支配排序,設計了多目標帝國競爭算法,求解多目標進場排序與調度模型的帕累托解。同時,引入覆蓋率指標、空間分布指標、超體積指標和平均理想距離等指標,衡量帕累托解集優劣。

4) 從兩個角度實施仿真驗證:一方面,利用基準案例(OR數據集)對比多目標帝國競爭算法(MOICA)與多目標遺傳算法(NSGA-II)及多目標模擬退火算法(MOSA)的解集質量與求解效率;另一方面,利用實際運行案例考察本文提出方法的有效性,并引入點融合系統保障優化結果的可行性。

5) 本文主要解決了單跑道進場排序與調度問題,未來的研究會擴展到多跑道、進離場、以及多機場終端空域的運行場景。

猜你喜歡
排序優化
排排序
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
排序不等式
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
恐怖排序
節日排序
刻舟求劍
兒童繪本(2018年5期)2018-04-12 16:45:32
主站蜘蛛池模板: 四虎综合网| 不卡视频国产| 99热这里只有精品在线播放| 亚洲国产清纯| 一级黄色片网| 黑色丝袜高跟国产在线91| 一级香蕉视频在线观看| 中文字幕人成乱码熟女免费| 天堂av高清一区二区三区| 欧美yw精品日本国产精品| 天堂网亚洲系列亚洲系列| 国产一区二区免费播放| 丁香六月综合网| 国产精品成人久久| 久久综合AV免费观看| 国产永久在线视频| 亚洲精品动漫| 免费av一区二区三区在线| 77777亚洲午夜久久多人| 欧美日韩国产在线播放| 亚洲天堂首页| 在线中文字幕日韩| 亚洲精品人成网线在线| 无码精品福利一区二区三区| www中文字幕在线观看| swag国产精品| 精品人妻无码区在线视频| 久久久成年黄色视频| 91福利一区二区三区| 亚洲成年人网| 国产黑丝一区| 国产色婷婷| 欧美日韩在线亚洲国产人| 国产永久无码观看在线| 亚洲中久无码永久在线观看软件| 999国产精品| 无码精油按摩潮喷在线播放| 在线观看国产精品第一区免费| 国产久草视频| 欧美中文字幕在线视频| 亚洲福利一区二区三区| 无码综合天天久久综合网| 精品久久久久久中文字幕女| 萌白酱国产一区二区| 国产美女免费网站| av尤物免费在线观看| 久久99久久无码毛片一区二区| 久久香蕉国产线看观看式| 国产爽妇精品| 国产啪在线| 精品欧美一区二区三区在线| 欧美国产综合视频| 青草视频在线观看国产| 国产成人a在线观看视频| 欧美日韩激情在线| 久久一本精品久久久ー99| 国产在线高清一级毛片| 国产情侣一区二区三区| 伊人成人在线| 草草影院国产第一页| 亚洲av无码牛牛影视在线二区| 精品国产美女福到在线直播| 欧美一级高清片欧美国产欧美| AV不卡国产在线观看| 99中文字幕亚洲一区二区| 99热这里只有精品在线播放| a级毛片毛片免费观看久潮| 人妻丰满熟妇啪啪| 在线观看国产黄色| 99精品视频在线观看免费播放| 久热中文字幕在线| 亚洲精品少妇熟女| 午夜啪啪网| 午夜a视频| 欧美日韩专区| 九色最新网址| 98超碰在线观看| 亚洲综合精品香蕉久久网| 亚洲欧美另类日本| 国产精品久久久久久久久| 国产18在线| 久久久亚洲色|