亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        柔性空域結(jié)構(gòu)下連續(xù)下降航跡多目標(biāo)優(yōu)化

        2021-03-26 09:52:32楊磊李文博劉芳子陳雨童趙征
        航空學(xué)報(bào) 2021年2期
        關(guān)鍵詞:優(yōu)化

        楊磊,李文博,劉芳子,3,陳雨童,趙征

        1. 南京航空航天大學(xué) 民航學(xué)院,南京 211106 2. 浙江大學(xué) 軟件學(xué)院,杭州 315100 3. 中國民用航空局 空中交通管理局 戰(zhàn)略發(fā)展部,北京 100020

        世界范圍內(nèi)空中交通系統(tǒng)正處于全面轉(zhuǎn)型升級階段,以應(yīng)對不斷增長的飛行需求、復(fù)雜龐大的體系結(jié)構(gòu)和多元多變的運(yùn)行環(huán)境,持續(xù)滿足各利益相關(guān)方對于飛行安全、運(yùn)行效率、成本效益和環(huán)境影響等航空運(yùn)營綜合性能期望[1]。基于航跡的運(yùn)行、基于性能的空中交通管理已成為全球航空運(yùn)輸系統(tǒng)變革的關(guān)鍵詞[2-4]。

        連續(xù)下降運(yùn)行(Continous Descending Operation, CDO)是提升終端區(qū)進(jìn)場交通流運(yùn)行效率和環(huán)保效益的關(guān)鍵技術(shù),最初由美國多個(gè)部門聯(lián)合組成的研究小組為減小航空噪音污染而設(shè)計(jì)[5]。之后,學(xué)者們開展了一系列的飛行實(shí)驗(yàn),證實(shí)了與傳統(tǒng)梯級下降程序相比,CDO通過使得航空器在降落前盡可能保持在較高高度,并盡可能使用最小發(fā)動(dòng)機(jī)推力和低阻力構(gòu)型下降,能夠有效減小燃油消耗、污染物排放及噪聲污染,具有顯著的經(jīng)濟(jì)和環(huán)保優(yōu)勢[6-9],引起了各國學(xué)者和業(yè)界的熱切關(guān)注。2010年,國際民用航空組織發(fā)布了Doc.9931《連續(xù)下降運(yùn)行手冊》,提供了CDO設(shè)計(jì)指南[10]。目前,CDO已經(jīng)成為歐洲單一天空計(jì)劃、美國下一代航空運(yùn)輸系統(tǒng)、國際民航組織航空系統(tǒng)組塊升級計(jì)劃,以及中國空管現(xiàn)代化發(fā)展戰(zhàn)略中的重點(diǎn)內(nèi)容之一。

        然而,研究指出高度剖面[11]、速度剖面[11]、下滑角[12]、氣象[13]、航班時(shí)刻[14]以及容量[15]等因素對于CDO運(yùn)行性能和環(huán)保性能具有較大影響,多機(jī)模式的CDO航跡優(yōu)化研究[16-17]成為新的研究熱點(diǎn)。學(xué)者們將基于時(shí)間管理(Time Based Management,TBM)和基于性能導(dǎo)航(Performance Based Navigation, PBN)的理念相融合[18],并通過探究航空器間隔[19],設(shè)置Metering Point[20],Required Time of Arrival[21-22]等方式來推動(dòng)CDO的實(shí)際應(yīng)用。此外,部分學(xué)者也研究了雷暴[23]、陣風(fēng)[24]等不確定因素影響下的CDO航跡實(shí)時(shí)更新,以及偏離優(yōu)化后航跡的實(shí)時(shí)調(diào)整策略等問題[25]。隨著研究和試點(diǎn)的不斷深入,傳統(tǒng)終端區(qū)標(biāo)準(zhǔn)進(jìn)離場程序較為剛性,難以支撐CDO最優(yōu)航跡規(guī)劃與實(shí)施的問題愈發(fā)凸顯。Alam等[26]提出一種新的過渡空域概念,以期通過構(gòu)建終端區(qū)內(nèi)的航跡調(diào)整緩沖區(qū),為動(dòng)態(tài)搜索最優(yōu)CDO航跡提供輔助。然而,該運(yùn)行概念與技術(shù)的靈活性有待提升,例如航空器一旦中止進(jìn)近就不能重新開始,并且在低高度區(qū)域需要不斷調(diào)整航向,增加了管制員和飛行員的工作難度和負(fù)荷[26-27]。此外,受目前地基管制自動(dòng)化系統(tǒng)以及空地航跡共享程度的制約,依賴管制員四維空間認(rèn)知的CDO模式僅在部分非繁忙機(jī)場或時(shí)段開展試點(diǎn),難以應(yīng)用于高密度終端區(qū)空域和復(fù)雜交通場景[17]。

        2010年,EUROCONTROL提出了點(diǎn)融合(Point Merge,PM)運(yùn)行概念和飛行程序,如圖1所示。PM主要由多個(gè)排序邊和一個(gè)融合點(diǎn)組成。航空器在排序邊上完成著陸排序,并在合適的時(shí)間調(diào)整航向和下降高度并直飛融合點(diǎn)。研究表明,PM扇形空域內(nèi)航班運(yùn)行具有較高的靈活性和航跡可預(yù)測性,有助于降低管制難度和工作負(fù)荷,提高飛行效率[28-31]。Errico等[32]發(fā)現(xiàn)基于點(diǎn)融合程序的CDO比傳統(tǒng)空域具有顯著的環(huán)保優(yōu)勢和運(yùn)行效益[33-35]。然而,由于點(diǎn)融合程序中排序邊的長度限制,其對延誤的吸收非常有限。高密度運(yùn)行時(shí),一旦達(dá)到點(diǎn)融合程序的飽和狀態(tài),需將航空器從點(diǎn)融合程序中向外引導(dǎo),對于空域內(nèi)航班運(yùn)行安全和效率具有顯著負(fù)向影響[36]。

        圖1 經(jīng)典點(diǎn)融合結(jié)構(gòu)圖Fig.1 Route structure of classic point merge

        因此,為進(jìn)一步提升CDO運(yùn)行能力和綜合效益,亟需打破傳統(tǒng)標(biāo)準(zhǔn)終端區(qū)儀表進(jìn)場航路(Standard Terminal Arrival Routes,STAR)剛性模式,設(shè)計(jì)更為柔性的空域結(jié)構(gòu)和靈活的運(yùn)行程序。本文以四維航跡運(yùn)行下星基導(dǎo)航模式為前提,考慮TBO下空中交通管制認(rèn)知復(fù)雜性,定義柔性進(jìn)場空域?yàn)橛梢幌盗锌臻g位置點(diǎn)構(gòu)成,且具備融合不同方向進(jìn)場流,提供靈活多樣飛行路徑選擇,以適應(yīng)動(dòng)態(tài)交通需求的大容量運(yùn)行空間,在保持管制員情景意識(shí)穩(wěn)定性的同時(shí),支撐CDO四維航跡安全高效運(yùn)行。鑒于此,本文基于已有研究成果,以點(diǎn)融合程序?yàn)榛A(chǔ),提出一種新型柔性進(jìn)場空域結(jié)構(gòu):倒皇冠形進(jìn)場空域(the Inverted Crown-Shaped Arrival Airspace,ICSAA),設(shè)計(jì)了航空器在新型空域中的運(yùn)行模式,建立了ICSAA下無沖突CDO航跡多目標(biāo)優(yōu)化方法,以期克服傳統(tǒng)點(diǎn)融合程序弊端,在一定程度上促進(jìn)CDO的應(yīng)用實(shí)踐。

        本文的主要工作包括:

        1) 構(gòu)建了一種新型的倒皇冠形進(jìn)場空域及其基本運(yùn)行模式。為打破傳統(tǒng)結(jié)構(gòu)性終端空域和程序化運(yùn)行方式對于CDO性能的限制,設(shè)計(jì)了一種倒皇冠形的柔性進(jìn)場空域,能夠融合各個(gè)方向的進(jìn)場交通流并生成初始CDO航跡,并通過規(guī)定兩種基礎(chǔ)運(yùn)行模式:下降高度(Mode V)和水平飛行(Mode H),保證航空器有序飛行,以簡化認(rèn)知難度,提高運(yùn)行效率和空域容量。

        2) 建立了無沖突CDO航跡多目標(biāo)優(yōu)化模型及元啟發(fā)式求解算法。基于ICSAA空域結(jié)構(gòu),研究面向(預(yù))戰(zhàn)術(shù)階段的CDO航跡優(yōu)化問題,即以油耗和飛行時(shí)間作為優(yōu)化目標(biāo),綜合考慮航空器動(dòng)力學(xué)性能限制,提出了無沖突CDO航跡多目標(biāo)規(guī)劃模型,設(shè)計(jì)了基于帶精英保留策略的非支配排序遺傳算法(NSGA-Ⅱ: Non-dominated Sort Genetic Algorithm-Ⅱ)的多目標(biāo)優(yōu)化高效求解算法,生成Pareto最優(yōu)解集。通過設(shè)計(jì)單架、低密度和高密度不同交通場景,驗(yàn)證了所提模型算法的合理性和有效性,論證ICSAA空域支撐高密度連續(xù)下降運(yùn)行潛力。

        1 新型進(jìn)場空域設(shè)計(jì)

        1.1 ICSAA結(jié)構(gòu)

        國際民航組織發(fā)布的Doc.9931《連續(xù)下降運(yùn)行手冊》指出,基于合適的空域、程序設(shè)計(jì)和空中交通管制許可,CDO能夠通過優(yōu)化航空器飛行剖面,建立低阻力構(gòu)型,降低發(fā)動(dòng)機(jī)推力,從而減少燃油消耗。本文基于當(dāng)前一般進(jìn)近空域水平和垂直范圍,轉(zhuǎn)變基于STAR的運(yùn)行理念,通過改善和擴(kuò)展傳統(tǒng)點(diǎn)融合程序,設(shè)計(jì)了一種新型柔性進(jìn)場空域結(jié)構(gòu):倒皇冠形進(jìn)場空域(ICSAA),如圖2所示。

        圖2 倒皇冠形進(jìn)場空域Fig.2 Structure of ICSAA

        ICSAA保留了傳統(tǒng)點(diǎn)融合程序中通過設(shè)置航路點(diǎn),并以“直飛”指令來引導(dǎo)飛行,從而減少飛行員與管制員通信量的優(yōu)勢。特別地,為解決傳統(tǒng)點(diǎn)融合排序邊對點(diǎn)融合程序運(yùn)行性能的影響,本文將排序空間擴(kuò)展為一個(gè)多層圓環(huán),融合點(diǎn)位于圓心。如圖2(a)所示,本示例中的ICSAA設(shè)置了5個(gè)不同高度的圓環(huán),高度間隔為300 m。最外層圓環(huán)高度為6 000 m,其上的航路點(diǎn)代表進(jìn)場點(diǎn);第5層圓環(huán)的高度為4 800 m,在其內(nèi)側(cè)還設(shè)置了相同高度的多層航路點(diǎn),其中每個(gè)航路點(diǎn)代表CDO的起始點(diǎn)。為保持進(jìn)近航空器之間的安全間隔,每個(gè)圓環(huán)的圓周上設(shè)置若干航路點(diǎn),相鄰航路點(diǎn)之間間隔為10°。各層圓環(huán)所處高度與半徑的比值等于tan 3°,即航空器在1~5層下降過程中的最大下滑角為3°;航空器到達(dá)第5層后,根據(jù)起始點(diǎn)選擇不同,保持3°~4°之間的下滑角飛向融合點(diǎn)。為確保航空器能夠截獲下滑道,融合點(diǎn)的高度設(shè)置為900 m。當(dāng)航空器飛到融合點(diǎn)之后,保持高度,截獲下滑道開始最后進(jìn)近。

        便于描述,本文對ICSAA內(nèi)的航路點(diǎn)逐一編號(hào)。除第5層圓環(huán)外,將每個(gè)圓環(huán)的航路點(diǎn)按照如圖2(b)逆時(shí)針方向標(biāo)記為1~36的序數(shù),并按照“(層數(shù),序數(shù))”格式進(jìn)行命名(例如,“(2,7)”表示第2層正北方向的點(diǎn))。航空器在到達(dá)第5層 圓環(huán)上的航路點(diǎn)之后,開始進(jìn)行連續(xù)下降進(jìn)近,并在融合點(diǎn)截獲下滑道。為確保截獲平順性,本文將最大切入角設(shè)置為60°,即圖2(b)中的包絡(luò)面的中心角為120°。因此,高度層4 800 m處的3個(gè)扇形環(huán)的最外側(cè)圓環(huán)編號(hào)為(5,1)~(5,13),中間一層航路點(diǎn)編號(hào)為(5,14)~(5,26),最內(nèi)側(cè)航路點(diǎn)編號(hào)為(5,27)~(5,39)。

        1.2 基本運(yùn)行模式

        ICSAA內(nèi),航空器飛行過程可以簡單分為2種:一是航空器在水平飛行階段時(shí),高度保持不變,按照規(guī)定的繞向依次飛往各個(gè)航路點(diǎn);二是航空器在進(jìn)場階段控制油門和飛機(jī)構(gòu)型等改變高度。因此,將ICSAA中的機(jī)動(dòng)飛行表示為2種模式:Mode H和Mode V[37],飛機(jī)在ICSAA內(nèi)的飛行活動(dòng)即可視為模式保持和切換。

        此外,為了保證航空器在ICSAA中運(yùn)行時(shí)CDO航跡的可預(yù)測性,航空器應(yīng)沿多層圓環(huán)空域的航路點(diǎn)逐層下降,且原則上不能穿過ICSAA內(nèi)部。為此,本文定義了Separation Z,如定義3所示。當(dāng)且僅當(dāng)航空器當(dāng)前所在航路點(diǎn)與下一層目標(biāo)航路點(diǎn)之間的水平間隔不大于給定值時(shí),航空器可切換為Mode V;否則航空器保持Mode H沿圓環(huán)圓周飛行??梢姡谕鈱优判颦h(huán)上,航空器具備Mode H和Mode V雙重運(yùn)行模式,直至到達(dá)第5層后保持Mode V執(zhí)行連續(xù)下降進(jìn)近。

        在多架航空器的場景中,情況將更加復(fù)雜。當(dāng)航空器以Mode V下降高度,飛向下一高度層的圓環(huán)時(shí),不僅需要滿足Separation Z約束,而且要保證在下降過程中與其他航空器之間沒有潛在飛行沖突;否則,應(yīng)以Mode H保持高度沿圓環(huán)飛行尋找滿足下降條件的目標(biāo)點(diǎn)和時(shí)機(jī)。當(dāng)然,航空器在沿環(huán)形圓周飛行時(shí),航空器主要通過微調(diào)速度與前機(jī)保持安全間隔。極端地,當(dāng)航空器在其性能范圍內(nèi)均無法通過改變高度和水平速度的方式保持與周邊航空器的安全間隔時(shí),等待程序啟用。

        綜上,與傳統(tǒng)點(diǎn)融合程序相比,ICSAA因其環(huán)形設(shè)計(jì),最外側(cè)的圓環(huán)能夠融合來自各個(gè)方向的進(jìn)場交通流,可在環(huán)形排序邊上容納更多航空器,降低排序難度,提升運(yùn)行效率,亦可一定程度上支撐終端區(qū)進(jìn)離場交通流分離、飛行限制區(qū)避讓和高空航路截彎取直。此外,當(dāng)航空器下降到第4層圓環(huán)上時(shí),可以選擇4 800 m高度(即第5層)上的任意航路點(diǎn),以不同下滑角開始執(zhí)行連續(xù)下降進(jìn)近,增加了航跡選擇靈活性?;诖耍挛膶⑨槍CSAA運(yùn)行特點(diǎn),權(quán)衡飛行時(shí)間和燃油消耗,研究高密度航班流下的多目標(biāo)無沖突CDO航跡規(guī)劃問題。

        2 問題描述

        本文將ICSAA結(jié)構(gòu)抽象為一個(gè)網(wǎng)絡(luò)G=(V,E)。其中,網(wǎng)絡(luò)中的節(jié)點(diǎn)(i,j)表示終端區(qū)航路點(diǎn)所屬層數(shù)與序數(shù)。根據(jù)ICSAA結(jié)構(gòu),Mode H、Mode V以及Separation Z約束,構(gòu)造鄰接矩陣Θ=e[(in,jn),(im,jm)]表征節(jié)點(diǎn)之間的有向連通性。若(in,jn)→(im,jm)連通時(shí),e[(in,jn),(im,jm)]=1;否則為+∞。因此,航空器α的路徑p可以表示為一系列連通節(jié)點(diǎn)的集合。此時(shí),ICSAA新型柔性網(wǎng)絡(luò)中連續(xù)下降進(jìn)近的航跡優(yōu)化問題可以表述為:根據(jù)未來一段時(shí)間窗內(nèi)的各進(jìn)場航班計(jì)劃進(jìn)場點(diǎn)和預(yù)計(jì)到達(dá)時(shí)間,根據(jù)航空器性能,權(quán)衡優(yōu)化油耗和飛行時(shí)間雙重目標(biāo),在滿足時(shí)效性要求前提下,在網(wǎng)絡(luò)G內(nèi)連通性約束下統(tǒng)籌規(guī)劃航班飛行路徑和過點(diǎn)時(shí)間。

        3 航空器四維軌跡建模

        為準(zhǔn)確評估ICSAA內(nèi)實(shí)施CDO的性能,本文建立四維航跡模型模擬航空器飛行。Glover和Lygeros從牛頓動(dòng)力學(xué)出發(fā),提出了點(diǎn)質(zhì)量模型(Point Mass Model, PMM)[38]。EUROCONTROL基于能量守恒定律提出了總能量模型(Total-Energy Model, TEM)[39]。在上述研究基礎(chǔ)上,針對點(diǎn)質(zhì)量模型進(jìn)行改進(jìn),以此來計(jì)算航空器連續(xù)下降運(yùn)行過程中各飛行參數(shù)。

        (1)

        式中:x、y、h為航空器在三維坐標(biāo)系中的位置;v為空速;τ為攻角;ψ為航向;W=(ω1,ω2,ω3)∈R3表示風(fēng)速;T為發(fā)動(dòng)機(jī)推力;D和L分別為空氣阻力和升力,表達(dá)式分別為

        (2)

        (3)

        其中:S為機(jī)翼參考面積;ρ為空氣密度;CD、CL分別為阻力系數(shù)和升力系數(shù),計(jì)算公式為

        CD=CD0+CD2×(CL)2

        (4)

        (5)

        式中:CD0和CD2為BADA中提供的計(jì)算系數(shù)。

        在模擬航空器運(yùn)動(dòng)的過程中,推力的獲取至關(guān)重要。本模型借助牛頓運(yùn)動(dòng)第2定律,將推力的求解轉(zhuǎn)化為加速度的求解,并利用BADA推薦的速度-高度剖面[11],采用如下公式變換,獲得加速度:

        (6)

        v(h)=-1.943×10-6·h2+0.034 99·h+73.34

        (7)

        為了使得模型更加接近航空器的真實(shí)運(yùn)行狀態(tài),本文使用了BADA中的大氣模型來模擬大氣特性(溫度、壓力、密度)隨高度的變化。式(2)~式(3)中的空氣密度ρ(kg/m3)是高度h的函數(shù),使用理想氣體定律由該高度的大氣壓強(qiáng)p和溫度Tem計(jì)算得到:

        圖3 速度-高度剖面曲線擬合Fig.3 Fitting curve of speed-height profile

        (8)

        式中:R為空氣的氣體常數(shù),其值為287.052 87 m2/(K·s2)。Tem0為平均海平面標(biāo)準(zhǔn)大氣溫度,其值為288.15 K;ΔTem為平均海平面上給定非標(biāo)準(zhǔn)大氣與國際標(biāo)準(zhǔn)大氣之間的溫度差,本文設(shè)置為0;β為低于對流層頂高度的ISA溫度梯度,其值為-0.006 5 K/m;p0為平均海平面標(biāo)準(zhǔn)氣壓。

        基于上述航空器性能模型式(1)~式(8),記為性能方程組M,當(dāng)航空器飛行起訖點(diǎn)(即ICSAA中的航路點(diǎn),分別為(is,js)和(ie,je))、飛行路徑p以及初始速度v確定后,航空器依據(jù)BADA推薦的速度-高度剖面所得到的加速度進(jìn)行速度控制,即可推定飛行過程中任意時(shí)刻t的航空器實(shí)時(shí)位置(x,y,h)。因此,航空器任意時(shí)刻位置可以表示為如下映射關(guān)系,可用于描述航空器四維軌跡,進(jìn)而用于軌跡優(yōu)化中航空器潛在飛行沖突判斷。

        H(t,(is,js),(ie,je),p,v,M)=(x,y,h)

        (9)

        4 CDO軌跡優(yōu)化建模

        4.1 決策變量

        為標(biāo)識(shí)航空器α是否經(jīng)過第i層第j個(gè)航路點(diǎn),引入0-1變量χ(α)(i,j),如式(10)所示:

        (10)

        進(jìn)一步,為標(biāo)識(shí)航空器飛越節(jié)點(diǎn)之間的連通性,引入變量表示航空器是否選擇該邊組成CDO路徑。設(shè)2個(gè)節(jié)點(diǎn)分別為(i1,j1)和(i2,j2),滿足e[(i1,j2),(i1,j2)]=1,則

        (11)

        4.2 目標(biāo)函數(shù)

        為了評估CDO的運(yùn)行成本和運(yùn)行效率,本文設(shè)計(jì)了2個(gè)優(yōu)化目標(biāo):航空器從進(jìn)場點(diǎn)到飛越融合點(diǎn)所花費(fèi)的燃油和飛行時(shí)間,如式(12)所示。需要說明的是,由于噪音影響與機(jī)場周圍區(qū)域的地理特征密切相關(guān),本文所設(shè)計(jì)的新型終端區(qū)暫未加入地形障礙物等因素,因此,在優(yōu)化目標(biāo)中暫不將噪音影響納入,可在后續(xù)研究特定機(jī)場應(yīng)用時(shí)加以考慮[40]。

        (12)

        式中:α表示不同的航空器;f(α)和t(α)表示在每架航空器在各條路徑上所花費(fèi)的油耗和時(shí)間;N(i,j)表示與(i,j)連通的航路點(diǎn)。

        在計(jì)算油耗時(shí),采用了BADA中的TSFC(Thrust Specific Fuel Consumption)模型。對于噴氣式發(fā)動(dòng)機(jī):

        (13)

        式中:Cf1、Cf2、Cfcr為BADA給定的燃油流率計(jì)算參數(shù);Time為個(gè)體適應(yīng)度;fc為航空器的燃油流率,是一個(gè)瞬時(shí)值;則油耗為燃料消耗率在下降時(shí)段內(nèi)的積分Fc,可近似為從進(jìn)場點(diǎn)到FAF整個(gè)航跡上每段相同時(shí)間區(qū)間內(nèi)油耗離散值之和,本文設(shè)置4 s。

        4.3 約束條件

        1) 唯一性約束

        航空器α的路徑必須確保唯一性,即除融合點(diǎn)以外,航空器經(jīng)過航跡中任意一點(diǎn)(ip,jp)后的下一個(gè)航路點(diǎn)有且僅有一個(gè)。

        (ip,jp)∈{(i,j)|χ(α)(i,j)=1}

        (14)

        2) 程序性約束

        為了確保航空器α必須經(jīng)過ICSAA中的每一個(gè)高度層逐層下降,因此,每一層的χ(α)(i,j)的累加和都要大于或等于1:

        (15)

        3) 連通性約束

        航空器α所選的路徑必須符合ICSAA運(yùn)行規(guī)則,即航空器所選路徑的前后相鄰的航路點(diǎn)的在鄰接矩陣中的值為1。

        e[(i,j),N(i,j)]=1

        (16)

        4) 加速度約束

        加速度是一個(gè)與航空器性能相關(guān)的值,影響著航空器安全運(yùn)行所需的推力,對于加速度的約束,也在間接約束了航空器運(yùn)行所需的推力??紤]到旅客舒適度,加速度不能過大。BADA手冊建議最大縱向加速度:

        (17)

        5) 最后進(jìn)近速度約束

        當(dāng)航空器開始著陸時(shí),速度應(yīng)盡可能使自己截獲下滑道。但有一個(gè)最小速度限制。最小著陸速度Vmin通過式(18)計(jì)算:

        (18)

        式中:常數(shù)1.3為BADA手冊中為所有航空器運(yùn)行推薦的數(shù)值;Vstall為航空器在參考質(zhì)量mref下的失速速度。

        6) 安全間隔約束

        C(t,α1,α2)=

        (19)

        綜上,新型空域下CDO航跡多目標(biāo)優(yōu)化模型如式(19)所示:

        (20)

        4.4 算法設(shè)計(jì)

        此類路徑分配問題的復(fù)雜度隨著航空器架次增長以幾何級數(shù)增加,無法在多項(xiàng)式時(shí)間內(nèi)得到和驗(yàn)證解,屬于NP-Hard問題,且解之間存在獨(dú)立性,故而采用智能優(yōu)化算法來求解[41]?;诰⒈A舨呗缘姆侵渑判蜻z傳算法(NSGA-II)是求解此類型多目標(biāo)問題的經(jīng)典算法,在求解質(zhì)量和收斂效率方面具有綜合優(yōu)勢[42-44]。NSGA-II算法執(zhí)行過程如圖4所示。

        圖4 NSGA-Ⅱ算法步驟Fig.4 NSGA-Ⅱ algorithm procedure

        針對本問題的NSGA-Ⅱ適應(yīng)性改進(jìn)如下:

        1) 編碼與解碼

        2) 初始種群

        初始種群作為啟發(fā)式算法的迭代起點(diǎn),是影響種群進(jìn)化結(jié)果和算法效率的重要因素。在單架航空器路徑優(yōu)化時(shí),隨機(jī)生成初始解,形成初始種群,進(jìn)行路徑優(yōu)化,獲得Pareto前沿解。而在進(jìn)行多架航空器路徑優(yōu)化時(shí),采用直接隨機(jī)生成初始解的方法過于低效,故而本文隨機(jī)選取單架航空器的Pareto前沿解來拼接組合成為多架航空器進(jìn)場路徑優(yōu)化的個(gè)體,繼而形成初始種群,以此提高求解效率。

        3) 非支配排序和擁擠距離排序

        產(chǎn)生初始種群后,采用非支配排序方法對種群內(nèi)的所有個(gè)體進(jìn)行快速分層,形成多個(gè)不同等級的Pareto前沿,并根據(jù)擁擠比較算子來實(shí)現(xiàn)不同等級的非支配解的擁擠距離排序,以保持和增加解的多樣性。由非支配等級和擁擠度決定的適應(yīng)度更強(qiáng)的個(gè)體被移動(dòng)到種群前沿,用于產(chǎn)生下一代。

        4) 遺傳策略

        采用二進(jìn)制錦標(biāo)賽算法和精英保留策略對父代種群的個(gè)體進(jìn)行選擇,選出精英個(gè)體進(jìn)行交叉和變異得到子代種群。個(gè)體之間在交叉時(shí),采用基于航班的交叉方式,即以一架航空器所對應(yīng)的決策變量為單位進(jìn)行交叉,交叉位置隨機(jī)選擇;而在變異時(shí),對于每個(gè)基因位隨機(jī)選擇采用基本位變異。不滿足約束條件的個(gè)體將被賦予一個(gè)較大的值,并在求解過程中被自動(dòng)拒絕。結(jié)合父代群體和子代群體,計(jì)算各染色體的適應(yīng)度值,使用精英保留策略保留更優(yōu)的染色體來產(chǎn)生新的父代群體。

        5) 終止條件

        通過設(shè)置最大迭代次數(shù)和檢驗(yàn)每一代種群的收斂結(jié)果兩種方式來終止算法。其中,當(dāng)種群中Pareto最優(yōu)前沿解的數(shù)量超過種群規(guī)模的80%以上時(shí),算法停止。

        5 實(shí)驗(yàn)設(shè)計(jì)

        5.1 參數(shù)設(shè)置

        實(shí)驗(yàn)使用Intel Core i7-8550U CPU 1.80 GHz 四核處理器,16 GB內(nèi)存,Windows 10 操作系統(tǒng)的筆記本電腦,采用MATLAB語言進(jìn)行編程求解。NSGA-II算法控制參數(shù)設(shè)計(jì)如下:種群數(shù)量為2 000,最大迭代次數(shù)為80,變異概率為0.7, 交叉概率為0.3。

        假設(shè)重力勢能高度等于地面高度,真空速等于地速,國際標(biāo)準(zhǔn)大氣壓,靜風(fēng)。選取民航業(yè)主流機(jī)型——空客A320作為實(shí)驗(yàn)機(jī)型。BADA 中A320的基礎(chǔ)參數(shù)如表1所示。

        表1 空客A320性能參數(shù)Table 1 Performance parameters of Airbus 320

        5.2 算法對比

        為驗(yàn)證NSGA-Ⅱ?qū)τ诒疚倪B續(xù)下降軌跡規(guī)劃問題的求解優(yōu)越性,選取了速度受限的多目標(biāo)粒子群算法(Speed-constrained Multi-objective Particle Swarm Optimization,SMPSO)和強(qiáng)Pareto支配算法(Improving the Strength Pareto Evolutionary Algorithm, SPEA2:)2種當(dāng)前性能較優(yōu)且應(yīng)用廣泛的多目標(biāo)優(yōu)化算法[45],以油耗和飛行時(shí)間為坐標(biāo)軸,建立二維目標(biāo)函數(shù)空間,采用Generational Distance (GD)、Inverted Generational Distance (IGD)、HyperVolume (HV)、Spacing和Maximum Pareto Front Error (MPFE)[45]5類評價(jià)指標(biāo),如式(21)所示,全面論證算法期望性能,即以用較少的計(jì)算資源得到覆蓋整個(gè)搜索空間、分布均勻且逼近Pareto前沿的非支配解集。其中,GD、IGD、Spacing以及MPFE是以距離為基礎(chǔ),其值越小越好;而HV衡量非支配解在二維目標(biāo)函數(shù)空間覆蓋的區(qū)域面積,其值越大越好。需要注意的是,在解算下列評價(jià)指標(biāo)時(shí),應(yīng)對數(shù)據(jù)進(jìn)行歸一化。

        (21)

        為對比NSGA-Ⅱ、SMPSO和SPEA2這3類算法在不同問題規(guī)模下的求解性能特征,設(shè)計(jì)了低中高3種不同密度的交通場景,即1/6/12架航空器同時(shí)進(jìn)場。為增強(qiáng)可對比性,3類算法具有相同的收斂條件,并且,NSGA-Ⅱ和SPEA2的種群數(shù)量均為2 000,最大迭代次數(shù)均為80,交叉概率均為0.3,變異概率均為0.7,其交叉和變異策略保持一致;SMPSO算法中,粒子數(shù)量和Archive Size為2 000,最大迭代次數(shù)為80,最優(yōu)求解性能所對應(yīng)的粒子位置迭代更新參數(shù)設(shè)定如下:初始速度影響因子w=0.2,個(gè)體最優(yōu)影響因子C1=0.2,全局最優(yōu)影響因子C2=0.5。在上述基準(zhǔn)參數(shù)下,3類算法的求解性能如表2所示。

        由結(jié)果可知,盡管3類算法性能均隨著問題規(guī)模的增大而降低,但NSGA-Ⅱ在各類性能指標(biāo)中始終優(yōu)于其他2種算法。此外,以單架航空器為例,3類算法收斂表現(xiàn)如圖5所示。同樣,NSGA-Ⅱ的收斂速度明顯優(yōu)于其他2種算法。因此,針對本文研究對象,NSGA-Ⅱ能夠更為高效和準(zhǔn)確地收斂到Pareto最優(yōu)解集附近。

        表2 3類算法性能對比Table 2 Performance comparison of three different algorithms

        基于NSGA-Ⅱ算法,下面將分別從單架、多架同時(shí)進(jìn)場和多架連續(xù)進(jìn)場視角,深入探討新型柔性空域下的連續(xù)下降航跡優(yōu)化問題。

        5.3 單架航空器航跡優(yōu)化分析

        單架航空器場景中,以第1層圓環(huán)5號(hào)點(diǎn)為進(jìn)場點(diǎn)為例開展優(yōu)化驗(yàn)證,求解時(shí)間約為20 s。每一代的Pareto前沿隨著迭代次數(shù)的增加而不斷收斂,如圖5所示。優(yōu)化后的部分可行解如圖6所示,其中紅色菱形中的解為Pareto最優(yōu)前端,其所對應(yīng)的CDO軌跡如圖7所示。進(jìn)一步,為外層所有進(jìn)場點(diǎn)尋求最優(yōu)航跡解集,運(yùn)行時(shí)間均小于30 s。表3展示了各個(gè)進(jìn)場點(diǎn)經(jīng)過優(yōu)化以后的Pareto最優(yōu)前端中油耗最少和飛行時(shí)間最短的解。其中,表中路徑序列包括4個(gè)數(shù)字,分別表示航空器從第2高度層下降至第5層所經(jīng)過的航路點(diǎn)編號(hào)。

        圖5 不同算法單架航空器Pareto前端收斂對比Fig.5 Convergence comparison of Pareto fronts evolving with generation in single aircraft scenario

        通常情況下,飛行時(shí)間越長,燃油消耗越大,但是實(shí)驗(yàn)結(jié)果表明油耗與飛行時(shí)間存在一定的權(quán)衡空間。鑒于此,本文首先對航空器下降過程中各飛行參數(shù)展開分析。如圖7所示,航跡1按照5-5-5-5的路線飛行,對應(yīng)Pareto最優(yōu)前沿解集中油耗最小解;航跡2按照5-5-5-31的路線飛行,對應(yīng)Pareto最優(yōu)前沿解集中飛行時(shí)間最短解。顯然,結(jié)果差異來自第5個(gè)高度層航路點(diǎn)選擇。如ICSAA結(jié)構(gòu)所示,從第5層的不同航路點(diǎn)下降,意味著不同的CDO航跡下滑角。圖8(a)表示航空器高度隨時(shí)間的變化。可以看出,航跡2的高度總是高于航跡1。由于實(shí)驗(yàn)中飛機(jī)的速度均符合BADA推薦的標(biāo)準(zhǔn)速度剖面,速度與高度成正比,因此,航跡2速度大于航跡1 (如圖8(b)所示),花費(fèi)下降時(shí)間更短。

        圖6 單架航空器優(yōu)化結(jié)果Fig.6 Trajectory optimization in single aircraft scenario

        圖7 Pareto最優(yōu)解航跡圖Fig.7 Trajectories of Pareto optimal solutions

        表3 各進(jìn)場點(diǎn)理想路徑Table 3 Optimal route for each entry point

        最底層的環(huán)形排序邊高度為4 800 m,經(jīng)過該高度層后,航空器開始進(jìn)行連續(xù)下降進(jìn)近。從圖8(a)中可以看出,航跡2開始實(shí)施CDO的時(shí)刻大約在250 s,而航跡1開始實(shí)施CDO的時(shí)刻大約在150 s。因此,航跡2在環(huán)形區(qū)域停留的時(shí)間更長。在環(huán)形區(qū)域停留期間需要提供額外的推力,故油耗較大。一旦航空器開始實(shí)施CDO時(shí),發(fā)動(dòng)機(jī)處于慢車或接近慢車狀態(tài),燃油消耗接近于零。因此,航跡2的燃油消耗更多,整個(gè)飛行過程中的推力和燃油流率隨時(shí)間的變化如圖8(c)和圖8 (d)所示。

        圖8 不同下降航跡飛行動(dòng)力中特征對比Fig.8 Comparison of flight dynamics for different descending trajectories

        5.4 多航空器航跡優(yōu)化分析

        維持高密度終端區(qū)內(nèi)安全高效運(yùn)行始終是CDO應(yīng)用面臨的重要挑戰(zhàn)之一。為了進(jìn)一步探究基于ICSAA的高密度終端區(qū)CDO性能,本文分別設(shè)計(jì)了同時(shí)進(jìn)場和持續(xù)進(jìn)場兩類交通場景,研究ICSAA的瞬時(shí)和持續(xù)運(yùn)行能力??紤]到CDO下的管制沖突認(rèn)知難度,為盡可能避免由于人為介入所引發(fā)的安全和效率降低問題,本文所有生成的最優(yōu)航跡均無沖突。

        5.4.1 場景1:多航空器同時(shí)進(jìn)場

        為探究基于ICSAA的CDO運(yùn)行瞬時(shí)性能,本文設(shè)置了航空器數(shù)量4/6/9/12遞增的同時(shí)進(jìn)場交通需求,其進(jìn)場點(diǎn)等間隔均勻分布于ICSAA最外層圓環(huán)。每個(gè)交通場景下,基于航空器飛越最外層進(jìn)場點(diǎn)時(shí)間和速度,利用本文所提出的優(yōu)化模型為每架航空器規(guī)劃最佳路徑。優(yōu)化求解時(shí)間在1~3 min之間,Pareto最優(yōu)前沿如圖9所示。

        總體上,隨著航空器數(shù)量遞增,Pareto最優(yōu)前沿向右上方移動(dòng),平均油耗和平均飛行時(shí)間逐步遞增,且呈現(xiàn)逐漸平緩的趨勢,表明隨著空域擁堵程度的加深,油耗與飛行時(shí)間的權(quán)衡關(guān)系具有相對明顯的收斂特征。

        圖9 多架航空器同時(shí)進(jìn)場Pareto前端Fig.9 Pareto optimal fronts for multi-aircraft scenarios in which aircraft simultaneously arrive at entry point

        以4架航空器為例,繪制了Pareto最優(yōu)前沿兩端的解對應(yīng)的軌跡,如圖10所示。對比圖10(a)與圖10(b)發(fā)現(xiàn),除了紫紅色航跡外,所有航跡均不相同,且差別主要仍在于第5層航路點(diǎn)的選擇,即下滑角不同,與單機(jī)運(yùn)行類似。在多機(jī)場景中,通過選擇不同下滑角有利于解決航空器之間的沖突,但對CDO航跡的油耗和飛行時(shí)間具有較為顯著的影響。

        圖10 優(yōu)化航跡對比Fig.10 Comparison of optimized trajectories

        為進(jìn)一步評估優(yōu)化效果,本文選取優(yōu)化后的Pareto前沿兩端分別與單架航空器理想運(yùn)行結(jié)果進(jìn)行對比,如表4所示。需要說明的是,為增加結(jié)果的可對比性,多航空器油耗最小解中的理想飛行時(shí)間是單架航空器獨(dú)立飛行最低油耗對應(yīng)的飛行時(shí)間;同理,多航空器飛行時(shí)間最小解中的理想油耗是單架航空器獨(dú)立飛行最小飛行時(shí)間對應(yīng)的油耗,均可從表3獲取。

        結(jié)果表明,如表4所示,隨著航空器數(shù)量增大,平均油耗和平均運(yùn)行時(shí)間一致性增長,油耗與飛行時(shí)間的增加比例整體亦呈增長趨勢,且油耗的增加趨勢更加明顯,如圖11所示。在不同擁擠程度下,取Pareto最優(yōu)前端中油耗最小值時(shí),最小平均油耗增加比例不超過20%,對應(yīng)的飛行時(shí)間增加比例均小于5%;取Pareto最優(yōu)前端中飛行時(shí)間最小值時(shí),最小平均飛行時(shí)間增加比例不超過9%,對應(yīng)的油耗增加量均小于13%。特別注意的是在部分場景下,油耗出現(xiàn)了負(fù)增加,其原因在于,本文選擇了單架航空器的Pareto最優(yōu)前沿兩端作為參考,當(dāng)航空器的飛行時(shí)間最小時(shí),油耗相對較大;而在多架航空器的場景下,由于航空器必須要避免沖突,一般不能采用單架航空器的Pareto最優(yōu)前沿中飛行時(shí)間最小解所對應(yīng)的路徑飛行,在交通密度較低時(shí)存在油耗更低的路徑解空間。然而,隨著航空器架次和沖突頻率的增加,航空器必須采取更多的機(jī)動(dòng)飛行以避免沖突,飛行時(shí)間和油耗呈現(xiàn)一致性增大,并均超過單架航空器Pareto最優(yōu)前端性能。

        綜上,在同時(shí)進(jìn)場交通場景下,本文所提出的

        表4 多架航空器同時(shí)進(jìn)場優(yōu)化結(jié)果對比

        圖11 相較理想解增加比例對比(同時(shí)進(jìn)場)Fig.11 Increased proportion of fuel consumption and trip time compared with ideal ones (for simultaneous arrivals)

        新型空域下CDO航跡優(yōu)化方法能夠權(quán)衡飛行時(shí)間與燃油消耗,為各航空器生成無沖突四維軌跡;與此同時(shí),前沿解的額外油耗和額外飛行時(shí)間均在可接受范圍內(nèi),具有應(yīng)對瞬態(tài)高密度交通需求的穩(wěn)定高效性能。

        5.4.2 場景2:多航空器持續(xù)進(jìn)場

        為探究基于ICSAA的CDO持續(xù)運(yùn)行性能,本文設(shè)計(jì)了進(jìn)場點(diǎn)固定,到達(dá)時(shí)間概率生成,數(shù)量不斷遞增的交通場景。設(shè)置了4個(gè)固定的進(jìn)場點(diǎn),編號(hào)分別設(shè)置為1、10、19和28。同一進(jìn)場點(diǎn)相鄰航空器之間的到達(dá)時(shí)間間隔服從泊松分布,數(shù)學(xué)期望為2 min。同時(shí),為了保證安全,相鄰兩架航空器之間仍應(yīng)保持最小安全間隔。

        由于航空器初始狀態(tài)與優(yōu)化效果緊密關(guān)聯(lián),因此,為評估航空器在極限狀態(tài)下的優(yōu)化效果,本文基于表3中油耗最低路徑為各航空器生成經(jīng)濟(jì)計(jì)劃航跡,采用Monte Carlo仿真(1 000 000次)提取沖突強(qiáng)度最大的交通場景進(jìn)行驗(yàn)證分析。其中,沖突強(qiáng)度計(jì)算采用離散時(shí)間(間隔4 s)沖突探測方法,統(tǒng)計(jì)全部進(jìn)場航空器在飛行過程中所遭遇的潛在沖突次數(shù)之和。

        航空器數(shù)量以4為單位從4架/20 min增加到16架/20 min,Monte Carlo仿真下的最大沖突強(qiáng)度如圖12所示,利用所提模型及算法開展航跡優(yōu)化,求解時(shí)間約為3~5 min,優(yōu)化結(jié)果如圖13所示??傮w上,隨著空域密度的增加,飛行效率與油耗之間的權(quán)衡空間先增大后減小。低密度場景下,盡管航空器潛在飛行沖突較少,但由于飛行速度采用最優(yōu)的速度-高度剖面,且航班數(shù)量較少,最優(yōu)航跡組合下可權(quán)衡空間相對較小。高密度場景下,航空器之間的運(yùn)行耦合程度高,空域態(tài)勢復(fù)雜,無沖突航跡解空間較小,最優(yōu)前沿緊縮。特別地,在中等密度條件下,如進(jìn)場率為12架次/20 min 時(shí),航跡優(yōu)化存在較為顯著的權(quán)衡區(qū)間,更能體現(xiàn)新型空域下CDO航跡規(guī)劃的彈性空間。

        圖12 最復(fù)雜交通場景沖突強(qiáng)度Fig.12 Collision intensity in simulated traffic scenarios with highest complex

        圖13 多架航空器連續(xù)進(jìn)場Pareto前端Fig.13 Pareto optimal fronts for trajectory optimization in continuous arrival scenarios

        同樣,取Pareto最優(yōu)前沿兩端與理想運(yùn)行結(jié)果進(jìn)行對比分析。隨著航空器數(shù)量的增加,優(yōu)化后的油耗和飛行時(shí)間增加比例呈現(xiàn)增長趨勢,如圖14所示。從最小油耗角度,不同密度的航空器運(yùn)行時(shí),最小油耗增加比例不超過30%,其所對應(yīng)的平均飛行時(shí)間不超過10%。從最小飛行時(shí)間角度,最小平均飛行時(shí)間增加比例不超過10%,其所對應(yīng)的油耗增加比例不超過15%,詳見表5。

        表5 多架航空器持續(xù)進(jìn)場優(yōu)化結(jié)果

        綜上,在持續(xù)進(jìn)場場景下,本文所提新型空域下的航跡優(yōu)化算法依然高效可靠,盡管隨著交通密度的增加,優(yōu)化后額外平均飛行時(shí)間和油耗比例持續(xù)增長,但從增長絕對量上依然處于高效可接受水平,與理想性能相比,額外平均油耗36~72 kg,平均飛行延誤約1~2 min,能夠有效支撐高密度終端區(qū)高性能運(yùn)行要求,有助于促進(jìn)TBO模式下空管和航空公司的協(xié)同性水平。

        圖14 相較理想解增加比例對比(持續(xù)進(jìn)場)Fig.14 Increased proportion of fuel consumption and trip time compared with ideal ones (for continuous arrivals)

        6 結(jié) 論

        1) 深入剖析了CDO運(yùn)行的限制因素,受傳統(tǒng)點(diǎn)融合技術(shù)的啟發(fā),突破現(xiàn)有的終端區(qū)空域結(jié)構(gòu)和運(yùn)行模式,構(gòu)建了支持全向進(jìn)場的一種倒皇冠形柔性空域,并初步定義了航空器在該空域內(nèi)的簡潔運(yùn)行方式。

        2) 以航空器飛行時(shí)間和油耗作為優(yōu)化目標(biāo),綜合考慮各類約束限制,建立CDO航跡多目標(biāo)優(yōu)化模型,改進(jìn)設(shè)計(jì)NSGA-Ⅱ高效求解算法,能夠支撐戰(zhàn)術(shù)階段無沖突四維航跡規(guī)劃。

        3) 對比驗(yàn)證了NSGA-Ⅱ算法優(yōu)越性,全面探討了所提模型算法在單架、多架同時(shí)進(jìn)場和多架持續(xù)進(jìn)場等各類場景下的優(yōu)化效果和特征,驗(yàn)證了ICSAA內(nèi)CDO航跡優(yōu)化應(yīng)對復(fù)雜交通態(tài)勢的有效性和高效性,對于支撐高密度終端區(qū)CDO運(yùn)行具有顯著潛力,有助于推動(dòng)TBO概念的應(yīng)用落地。

        4) 本論文是對新型空域下CDO航跡優(yōu)化的理論探索,后續(xù)將采用中國機(jī)場及空域?qū)嶋H數(shù)據(jù),進(jìn)一步研究ICSAA在不同空域限制下的適應(yīng)性變構(gòu)技術(shù),將噪音影響納入航跡多目標(biāo)優(yōu)化中。此外,如何建立基于動(dòng)態(tài)空地協(xié)商的最優(yōu)航跡多準(zhǔn)則決選,也是下一步所開展的重要方向,為TBO下航跡協(xié)同和自主協(xié)商提供輔助決策。

        [21] PAWELEK A, LICHOTA P, DALMAU R, et al. Fuel-efficient trajectories traffic synchronization[J]. Journal of Aircraft, 2019, 56(2): 481-492.

        猜你喜歡
        優(yōu)化
        超限高層建筑結(jié)構(gòu)設(shè)計(jì)與優(yōu)化思考
        PEMFC流道的多目標(biāo)優(yōu)化
        能源工程(2022年1期)2022-03-29 01:06:28
        民用建筑防煙排煙設(shè)計(jì)優(yōu)化探討
        關(guān)于優(yōu)化消防安全告知承諾的一些思考
        一道優(yōu)化題的幾何解法
        由“形”啟“數(shù)”優(yōu)化運(yùn)算——以2021年解析幾何高考題為例
        圍繞“地、業(yè)、人”優(yōu)化產(chǎn)業(yè)扶貧
        事業(yè)單位中固定資產(chǎn)會(huì)計(jì)處理的優(yōu)化
        4K HDR性能大幅度優(yōu)化 JVC DLA-X8 18 BC
        幾種常見的負(fù)載均衡算法的優(yōu)化
        電子制作(2017年20期)2017-04-26 06:57:45
        亚洲写真成人午夜亚洲美女| 久久天天躁狠狠躁夜夜av| 国模欢欢炮交啪啪150 | 亚洲另类欧美综合久久图片区 | 午夜理论片yy6080私人影院| 女同久久精品国产99国产精品 | 成年男女免费视频网站| 亚洲欧美一区二区三区国产精| 最新国产主播一区二区| 日韩一级137片内射视频播放| 国产高清在线一区二区不卡| 国产午夜免费高清久久影院| 国产成人午夜精华液| 亚洲免费观看| 日本熟妇高潮爽视频在线观看| 国产一区av男人天堂| 欧美三级不卡在线观看| 性一交一乱一伦| 97人妻无码免费专区| 亚洲一区二区日韩精品| 肥老熟女性强欲五十路| 99精品国产一区二区| 一本久道久久综合五月丁香| 一区二区三区观看在线视频| 久久无码潮喷a片无码高潮| 亚洲av永久无码精品网址| 免费人成视频在线观看网站| 欧美性xxx久久| 日韩字幕无线乱码免费| 无码精品一区二区三区在线| 亚洲女初尝黑人巨高清| 亚洲熟妇色xxxxx欧美老妇| 一区二区三区视频免费观看在线| 久久久精品国产免大香伊| 人妻无码久久一区二区三区免费| 国产精品女丝袜白丝袜| 亚洲综合在线观看一区二区三区| 久久不见久久见中文字幕免费| 大陆一级毛片免费播放| 麻豆成年人视频在线观看| 在线观看免费无码专区|