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

基于形態(tài)特征的終端區(qū)進場中心航跡識別方法

2023-10-09 01:57:00王超李昊昱陳含露
科學(xué)技術(shù)與工程 2023年26期

王超, 李昊昱, 陳含露

(中國民航大學(xué)空中交通管理學(xué)院, 天津 300300)

空中交通航跡數(shù)據(jù)反映了真實的空中交通情況,記錄了真實的飛行路徑、高度航向變化等細(xì)節(jié),是管制員對大量航空器調(diào)度、指揮和引導(dǎo)的結(jié)果,空中交通軌跡的形態(tài)特征和空間分布特性是空中交通復(fù)雜性重要的外在表現(xiàn)形式。中心航跡是航空器飛行過程中具有代表性的軌跡,它描述了航空器群相似的運動趨勢。在復(fù)雜的終端區(qū)運行環(huán)境下通過對軌跡進行相似性度量,聚類并提取中心航跡,對表征交通流特征、挖掘航空器飛行模式和檢測異常軌跡[1],評估實際飛行軌跡與飛行程序的一致性[2],研究客觀表征軌跡運動混亂程度的復(fù)雜性模型,量化分析空中交通系統(tǒng)性能等具有重要意義。

當(dāng)前已有諸多學(xué)者開展了軌跡聚類及相似性度量的相關(guān)研究。王超等[3]考慮飛行速度差異帶來的航跡點相似度的偏差,將雷達航跡點逆向比對,提出基于軌跡聚類的平均軌跡構(gòu)造方法。熊偉等[4]通過相異度矩陣對數(shù)據(jù)進行預(yù)聚類處理,使用高斯混合EM算法對軌跡進行聚類,提高了聚類的準(zhǔn)確性。楊璐等[5]建立了多特征融合的軌跡相似度模型,通過自適應(yīng)譜聚類算法對軌跡進行分析,以各簇軌跡間距離為指標(biāo)提取中心航跡。Besse等[6]對基于動態(tài)時間扭曲和基于形狀的兩種相似度度量方法的聚類結(jié)果進行了定量分析,發(fā)現(xiàn)基于形狀的距離更適用于軌跡聚類,聚類主要考慮軌跡的空間緊湊性,而時間和速度等非空間特征的影響并非十分重要。Wang等[7]提出了一種考慮軌跡形狀的相似性度量方法,但沒有考慮到軌跡的連續(xù)性。Li[8]考慮了豪斯多夫距離在計算軌跡間距時容易受到航跡點跨度大、航跡點缺失等問題的影響,利用改進的DBSCAN算法進行軌跡段聚類,通過提取質(zhì)心向量得到中心航跡。

現(xiàn)有軌跡聚類方法未考慮到軌跡整體的形狀、軌跡段方向、速度、連續(xù)性等特征,同時易受到噪聲點和軌跡長度的影響[9],無法在大量航跡數(shù)據(jù)中直觀地挖掘出軌跡的整體運動特征并精確識別出中心航跡。因此,現(xiàn)提出一種基于單向距離(one way distance, OWD)和密度峰值聚類(density-peak clustering, DPC)的中心航跡提取方法。根據(jù)軌跡空間形態(tài)特征,采用單向距離相似度度量方法,同時考慮軌跡的位置屬性和航向?qū)傩远x多特征軌跡相似度模型。通過密度峰值聚類進行分析,提取聚類結(jié)果中的每一簇中具有最高密度的真實軌跡作為中心航跡。最后通過實例分析,驗證所提模型和方法的合理性和有效性。

1 軌跡相似度

軌跡相似性度量是軌跡聚類的基礎(chǔ)[10]。軌跡相似度是指一對軌跡之間的相似程度,包括時空位置、形狀趨勢、運動特征等,通過軌跡相似度可以衡量運動的整體相似性。基于幾何形狀計算相似度的主要有豪斯多夫距離和弗雷歇距離[11-12],它們都是基于點的相似度計算方法,易受到噪聲較大的影響。現(xiàn)有的測量軌跡間相似性的方法多具有以下局限性:對噪聲和誤差非常敏感,魯棒性差;多考慮時間和空間因素,忽略了軌跡形狀、軌跡段方向等特征;在匹配軌跡時需要單調(diào)連續(xù);對采樣率和采樣個數(shù)要求高,前期需要做大量的數(shù)據(jù)預(yù)處理工作[9]。本文研究在形狀相似度度量方法單向距離的基礎(chǔ)上進行相似度函數(shù)的構(gòu)造,提出一種基于單向距離的多特征軌跡相似度模型。單向距離采用點到線段的距離而不是點到點之間的距離,分段對稱計算軌跡之間的相似性,將軌跡作為一個整體考慮,包括軌跡的形狀和物理距離,能夠?qū)Σ煌L度的路徑間距離進行歸一化,對軌跡中的噪聲點具有較強抗干擾能力,降低了對采樣率的要求,且無需軌跡單調(diào)連續(xù)。

1.1 單向距離

使用網(wǎng)格對軌跡離散化處理,將整個區(qū)域分為相同大小的單元網(wǎng)格,每個單元網(wǎng)格根據(jù)其在x和y坐標(biāo)中的位置進行標(biāo)記。假設(shè)兩個單元網(wǎng)格g1(x1,y1)和g2(x2,y2),則兩個單元網(wǎng)格之間的距離為

(1)

網(wǎng)格軌跡Tg=(g1,g2,…,gn)由單元網(wǎng)格序列組成,相鄰單元網(wǎng)格的距離為1,n為Tg的長度,記為|Tg|=n。

對于一個單獨的單元網(wǎng)格g和一條軌跡Tg,如果網(wǎng)格軌跡Tg中的網(wǎng)格g′∈T到網(wǎng)格g的距離D(g′,g)小于網(wǎng)格軌跡Tg中任意其他網(wǎng)格g″到網(wǎng)格g的距離D(g″,g),則將單元網(wǎng)格g′定義為網(wǎng)格軌跡Tg上的局部最小點。單元網(wǎng)格g到局部最小點的距離D(g,g′)即為其到軌跡的最短距離。

(2)

假設(shè)航跡數(shù)據(jù)集合中任意兩條軌跡分別為Ti和Tj。從軌跡Ti到軌跡Tj的單向距離Dowd(Tgi,Tgj)定義為網(wǎng)格軌跡Tgi的局部最小點p到網(wǎng)格軌跡Tgj的距離除以Tgi的長度,該距離是從網(wǎng)格軌跡Tgi到網(wǎng)格軌跡Tgj的最短距離的積分,為有向距離。

(3)

兩條軌跡Ti和Tj之間的單向距離Ds(i,j)是它們的單向距離之和的平均值,即得到的Ds(i,j)是對稱的。

(4)

1.2 多特征軌跡相似度模型

在終端空域中,航空器一般都是按照規(guī)定的標(biāo)準(zhǔn)飛行程序從不同的方向進場,僅使用單向距離的度量方法并不完全適合于處理空管航跡數(shù)據(jù),因此在使用單向距離計算軌跡相似度的基礎(chǔ)上,考慮到兩條軌跡之間的位置屬性和航向?qū)傩?定義軌跡之間的位置特征距離和航向特征距離。

航跡數(shù)據(jù)是一個多維序列,由多維數(shù)據(jù)點組成每個航跡點包括經(jīng)度、緯度、地速、高度和時間等信息[13]。假設(shè)航跡數(shù)據(jù)集合中第i條和第j條軌跡分別表示為

(5)

(6)

軌跡Ti與Tj的所有航跡點的歐式距離d={d1,d2,…,dk},定義軌跡Ti到軌跡Tj的位置特征距離為

(7)

定義軌跡Ti到軌跡Tj的航向特征距離為

Dθ(i,j)=|θi-θj|

(8)

式(8)中:θi和θj分別為軌跡Ti和Tj靠近終端區(qū)邊界處航跡點的航向平均值。

定義軌跡Ti到另一條軌跡Tj的多特征軌跡相似度模型為

D(i,j)=ωsDs(i,j)+ωdDd(i,j)+ωθDθ(i,j)

(9)

式(9)中:Ds(i,j)為軌跡Ti與Tj之間的單向距離;Dd(i,j)為軌跡Ti與Tj之間的位置距離;Dθ(i,j)為軌跡Ti和Tj之間的航向距離;ωs為單向距離的權(quán)重因子;ωd為位置距離的權(quán)重因子;ωθ為航向距離的權(quán)重因子。權(quán)重因子的取值取決于多特征距離的應(yīng)用場景,滿足ωs≥0,ωd≥0,ωθ≥0且ωs+ωd+ωθ=1。

(10)

軌跡之間的各特征距離進行歸一化處理,Xi為轉(zhuǎn)換前的數(shù)值,XN表示轉(zhuǎn)換后的數(shù)值,Xmax和Xmin分別為樣本最大值和最小值,nT表示軌跡總條數(shù)。經(jīng)過多特征軌跡相似度模型計算之后,得到對稱的距離矩陣DT,作為之后聚類分析的輸入。

(11)

2 密度峰值聚類分析及中心航跡提取

基于所定義的軌跡之間的相似度表達,通過無監(jiān)督聚類算法將軌跡按照空間相似度進行劃分。目前廣泛使用的聚類算法包括k均值算法、分層聚類、密度算法等[14]。線段簇通常是任意形狀的,基于密度的軌跡聚類方法可以發(fā)現(xiàn)任意形狀的聚類并且過濾噪聲,是最適合用于線段的聚類方法。密度峰值聚類(density-peak clustering, DPC)[15]基于兩個假設(shè):簇中心的局部密度大于圍繞它的其他數(shù)據(jù)點的局部密度;不同簇中心間的距離較遠(yuǎn)。基于多特征相似度模型的密度峰值聚類算法具體步驟及流程如下。

步驟 1對給定的航跡集,首先計算相似性距離矩陣DT,隨機選取k個初始對象(T1,T2,…,TK)作為初始聚類中心點。

步驟 2計算各個初始聚類中心的局部密度值ρi,并將各中心的密度從大到小排列。各軌跡的局部密度是以該軌跡為中心,截止距離dc為半徑的范圍內(nèi)的軌跡的數(shù)量,使用高斯核函數(shù)計算局部密度ρi,計算公式為

dc=dround[τ·mT]

(12)

(13)

式(13)中:D(i,j)為軌跡Ti和Tj的多特征相似度距離;dc為截止距離閾值參數(shù),將距離矩陣的距離dround按升序排列即為d1≤d2≤…≤dmT;mT為距離矩陣中距離的個數(shù);τ為百分比;[]為取整函數(shù)在以dc為半徑的范圍內(nèi)進行搜索計算。

步驟 3根據(jù)密度值計算初始聚類中心Ti與其最近的較高密度聚類中心Tj之間的上向距離δi,即

(14)

步驟 4將各樣本劃分到局部密度較大且與各聚類中心最近的類簇中,得到新的聚類中心,重復(fù)執(zhí)行步驟2和步驟3,至聚類中心不再變化。

步驟 5算法終止,輸出最終的聚類中心點和類簇。

從各類軌跡簇中計算并提取一條能夠代表整體運動方向的中心軌跡。計算得到局部密度最高和距離更高密度中心較遠(yuǎn)的軌跡,密度峰值聚類算法輸出的各聚類中心即為需要提取的中心航跡。

3 實例分析

3.1 實驗環(huán)境與數(shù)據(jù)預(yù)處理

以雙流國際機場終端區(qū)進場ADS-B航跡數(shù)據(jù)為例,使用MATLAB作為數(shù)據(jù)處理和實驗分析軟件,提取航空器運行的經(jīng)度、緯度、高度、速度、時間等航跡點信息,以及機型、起降機場等航班信息,對航跡點進行航跡整合、裁剪和清洗操作。將同航班號的航跡點數(shù)據(jù)按照時間順序構(gòu)成一條軌跡序列,刪除時間與經(jīng)緯度同時重復(fù)的相鄰航跡點。

空管監(jiān)視系統(tǒng)獲取的航跡數(shù)據(jù)盡管經(jīng)過相關(guān)、降噪、降重等處理,但是仍然包含大量冗雜無關(guān)的航跡信息。大部分航跡數(shù)據(jù)包含航路飛行階段,甚至是從起飛到落地的全部飛行過程。首先對軌跡數(shù)據(jù)進行裁剪,終端空域的地理范圍是一個多邊形,剔除那些超出終端空域范圍外的航跡,保留在多邊形上和多邊形內(nèi)部的航跡。在終端空域內(nèi)包含大量飛越航空器的航跡,也存在通航飛行、訓(xùn)練飛行的空域。在通航飛機裝備了ADS-B發(fā)射機的情況下,監(jiān)視系統(tǒng)同樣會記錄其軌跡。可以依據(jù)軌跡的起點、終點,高度變化和飛行范圍,通過判斷航空器是否在使用本場跑道落地來剔除這些航跡。

由于進離場航班存在完全不同的飛行模式,需對進離場飛行軌跡進行區(qū)分。依據(jù)軌跡對應(yīng)的航班飛行計劃數(shù)據(jù)提供的航班進離場時間、起落機場等信息,識別起飛機場與目的地機場,將軌跡數(shù)據(jù)劃分為進場與離場。對于經(jīng)停航班,由于前期已根據(jù)高度范圍將軌跡數(shù)據(jù)中地面滑行段航跡點剔除,計算軌跡中相鄰航跡點時間差,將軌跡劃分為進離場兩段軌跡。通過對軌跡數(shù)據(jù)的預(yù)處理,僅保留在終端空域范圍內(nèi)的進場商業(yè)航空器飛行過程中的航跡數(shù)據(jù)。

聚類時使用直角坐標(biāo)系進行計算,為提高數(shù)據(jù)處理與計算速度,選取終端區(qū)范圍內(nèi)航跡進行分析。雙流國際機場共有兩條跑道02L/20R以及02R/20L,如圖1所示,雙流國際機場主要使用IGNAK、EKOKA、AKDIK、LADUP、MEXAD方向進場程序進場,使用跑道02R的進場程序主要有8條,使用跑道02L的進場程序主要有7條,使用20L跑道進場程序較少,一般只因施工或天氣因素才使用。

圖1 雙流機場主要進場飛行程序示意圖

3.2 進場交通流識別及中心航跡提取

對航跡數(shù)據(jù)預(yù)處理后,基于多特征軌跡相似度模型和密度峰值聚類算法,以對雙流機場終端區(qū)單日414條進場飛行軌跡為研究對象,首先由式(1)~式(4)計算進場飛行軌跡的形狀距離,由式(6)~式(7)計算進場飛行軌跡的位置距離,之后通過式(8)計算進場飛行軌跡的航向距離。依據(jù)實際數(shù)據(jù)經(jīng)過多次試驗得到權(quán)重因子。通過式(9)對各個距離進行加權(quán),利用式(10)對所有距離進行歸一化,得到式(11)軌跡之間的所有軌跡相似性距離矩陣DT。DPC算法的參數(shù)截斷距離dc的值需要通過設(shè)置一個百分比τ來確定,τ取值范圍一般為1%~2%,確定截斷距離前,需要先根據(jù)式(12)對DT所有軌跡間的相似度距離D(i,j)做升序排列。本算例中過濾掉前1.4%的值,取此時的最小值為截斷距離。本算例中dc,計算結(jié)果為0.012 3。根據(jù)式(13)計算該軌跡的局部密度參數(shù)ρi。根據(jù)式(14)計算數(shù)據(jù)軌跡Ti到具有更高密度的軌跡的最近距離δi。當(dāng)局部密度和相對距離的值都比較大時,此時對應(yīng)的軌跡則為初始聚類中心即密度峰值點。

通過密度峰值聚類算法不斷迭代從而將軌跡劃分成不同類型的軌跡簇,進場交通流軌跡識別結(jié)果如圖2所示。該方法可以識別出符合標(biāo)準(zhǔn)進場程序的主要交通流和不同于進場程序的其他主要進場路徑,能較全面地反映出航空器進場軌跡的主要特征,從而掌握終端區(qū)航空器的飛行模式,分析軌跡簇的分布規(guī)律。實際運行中軌跡的時空分布十分復(fù)雜,在終端區(qū)內(nèi),航空器還存在著一些異常軌跡。產(chǎn)生異常飛行狀態(tài)一般由于管制員工作負(fù)荷大、管制水平低、工作經(jīng)驗不足,飛行員未執(zhí)行標(biāo)準(zhǔn)操作,空域飛行流量大和惡劣天氣等原因造成。將形態(tài)差異度過大、航向大幅改變、偏離標(biāo)準(zhǔn)程序、與大部分軌跡運動模式相異且數(shù)量較少的軌跡作為異常軌跡。使用密度峰值聚類算法得到單日終端區(qū)進場飛行軌跡各簇形態(tài)如圖3所示,識別出的8個類簇的軌跡彼此在形狀和空間分布上有著顯著的差異,而同一簇內(nèi)的軌跡有著較高的相似度和形態(tài)特征。統(tǒng)計得到不同飛行模式下的飛機數(shù)量,一定程度上反映出了終端區(qū)內(nèi)各交通模式的頻繁程度以及總體運行情況。

圖2 終端區(qū)進場飛行軌跡聚類結(jié)果

圖3 各簇進場飛行軌跡形態(tài)示意圖

根據(jù)實際情況對軌跡進行分類時,既需要保持簇內(nèi)軌跡的相似性,又要體現(xiàn)不同類簇軌跡的分布特征和差異性。類簇數(shù)目k對聚類的合理性具有重要的影響,當(dāng)類簇數(shù)目k的取值較小時,基本可以表達有特點的軌跡路線,但仍有部分較為典型的軌跡無法被表達。根據(jù)實際當(dāng)類簇數(shù)目k的取值較大時,可區(qū)分同一個進場且形態(tài)多樣的各類軌跡,但有些實質(zhì)為管制員引導(dǎo)航空器飛三邊等指揮路線的變形,應(yīng)視為同一類軌跡。在本算例中,當(dāng)k=8時效果最好,可以較為科學(xué)全面地反映出實際運行中航空器各類進場軌跡的主要特征。

對軌跡數(shù)據(jù)聚類分析得到軌跡簇后,為進一步挖掘各軌跡簇的整體運動規(guī)律,用一條或多條代表性的軌跡來描述海量軌跡的空間整體移動特征,這個過程就是中心航跡的提取。其他軌跡與這些中心航跡具有相同或相似的行徑路線。因此,通過計算聚類簇中的每條軌跡到簇中其他軌跡的相似度距離之和,將距離較遠(yuǎn)相似度較高且具有最高密度的一條軌跡作為中心航跡,在密度峰值聚類算法中,其聚類中心是具有局部最高密度且距離更高密度中心較遠(yuǎn)的軌跡,因此選取各簇軌跡的聚類中心作為中心航跡,如圖4所示,各軌跡簇聚類中心vi的上向距離和局部密度如表1所示。

表1 各軌跡簇聚類中心實驗參數(shù)表

圖4 單日終端區(qū)進場飛行軌跡各簇中心航跡圖

對航跡運行數(shù)據(jù)進行統(tǒng)計,實際數(shù)據(jù)表明,管制員只使用了部分偏好程序,雙流機場切入02跑道航向道的導(dǎo)航點位于距離跑道入口73.2 km處,然而在實際運行中航空器一般在距離跑道入口40 km處切入跑道航向道,這是出于安全和效率的考慮,在保證航空器安全間隔的條件下,內(nèi)移切入跑道航向道的導(dǎo)航點。從IGNAK方向進場的軌跡簇空間形態(tài)和位置特征與標(biāo)準(zhǔn)飛行程序差異顯著,該類軌跡整體直線段較多,可能是由于交通流量控制,空管運行部門直接指揮引導(dǎo)航空器直飛最后進近定位點、引導(dǎo)至某點再直飛五邊、飛三轉(zhuǎn)彎等距圓等實施機動飛行的指令造成的,可以一定程度上提高航空器進場的效率。

3.3 聚類效果評價

SD是基于密度的指標(biāo)[16],SD由簇間密度和簇內(nèi)方差組成,通過對比類內(nèi)的緊密性和類間的密度來評估聚類的有效性。SD的取值越小表明聚類效果越好,計算公式為

SD=SC+ρW

(15)

(16)

(17)

式中:ρW為類間的密度;SC為簇間的平均分離度;c為類簇數(shù)目;ci為第i類;ρ(vi)為第i類的聚類中心vi的密度;ρ(vj)為第j類的聚類中心vj的密度;ρ(uij)為第i類和第j類之間的中心uij的密度;uij為聚類中心vi和vj兩點連線的中點;‖δ(S)‖為整個樣本數(shù)據(jù)S的方差;‖δ(vi)‖為第i類樣本數(shù)據(jù)的方差。

輪廓系數(shù)(silhouette coefficient, SC)[5]通過內(nèi)聚度和分離度來評估聚類算法的性能,適用于實際類別信息未知的情況。假設(shè)某一軌跡和它位于的類簇Cc空間中的其他軌跡的相似度為ai,軌跡和其他類簇Cp中的軌跡的相似度為bi,則第i條軌跡的輪廓系數(shù)si可表示為

(18)

(19)

(20)

(21)

假設(shè)樣本在類簇C中;n為樣本總數(shù),nc為當(dāng)前類簇樣本數(shù),np為其他類簇樣本數(shù);輪廓系數(shù)越接近于1,聚類效果越理想。

對單日進場飛行軌跡聚類分析,將得到的中心航跡與軌跡簇通過SD指標(biāo)和輪廓系數(shù)指標(biāo)進行聚類效果評價,并與層次聚類算法進行對比,不同參數(shù)下的輪廓系數(shù)SC和SD指標(biāo)值如表2所示,可以看出,使用密度峰值聚類算法對雙流機場進場軌跡進行聚類分析時分為8個類簇時效果最好,且密度峰值聚類算法模型的評估結(jié)果要優(yōu)于層次聚類算法。

4 結(jié)論

(1)基于ADS-B航跡數(shù)據(jù),使用單向距離度量軌跡之間的相似性,根據(jù)軌跡的形狀、航向、位置特征構(gòu)建了相似性矩陣,利用密度峰值聚類算法對軌跡進行劃分并提取聚類質(zhì)心作為中心航跡。結(jié)果表明,模型能有效且直觀地識別出復(fù)雜環(huán)境下不同空間形態(tài)的軌跡簇。

(2)該方法為終端區(qū)航空器飛行模式挖掘,異常軌跡識別和終端空域空中交通時空復(fù)雜性及形成機理等相關(guān)研究提供了技術(shù)支持。接下來將圍繞評估實際運行軌跡與飛行程序的一致性,精確識別異常軌跡和探究廣義空中交通復(fù)雜性模型等做進一步的研究。

主站蜘蛛池模板: 在线播放国产99re| 91视频首页| 日韩精品少妇无码受不了| 国产黄色爱视频| 国内嫩模私拍精品视频| 动漫精品啪啪一区二区三区| 激情综合婷婷丁香五月尤物| av一区二区无码在线| 91九色视频网| 最新日本中文字幕| 久久国产黑丝袜视频| 九色综合视频网| 中字无码精油按摩中出视频| 欧美一级高清免费a| 18禁影院亚洲专区| 九色在线观看视频| 久久亚洲国产最新网站| 天天躁夜夜躁狠狠躁躁88| a亚洲天堂| 日韩毛片基地| 日韩精品一区二区三区免费| 国产女主播一区| 国产性生交xxxxx免费| 国内精品伊人久久久久7777人| 福利在线一区| 天堂网国产| 香蕉在线视频网站| 一本一道波多野结衣一区二区 | 精品综合久久久久久97| 精品视频一区二区观看| 69国产精品视频免费| 精品伊人久久大香线蕉网站| 91久久国产综合精品女同我| 第一页亚洲| 日韩高清成人| 999精品视频在线| 996免费视频国产在线播放| 99爱在线| 制服丝袜一区| 久久婷婷六月| 666精品国产精品亚洲| 亚洲国产综合精品一区| 色吊丝av中文字幕| 97亚洲色综久久精品| 午夜福利视频一区| 国产波多野结衣中文在线播放| 72种姿势欧美久久久大黄蕉| 92午夜福利影院一区二区三区| 亚洲国产精品一区二区第一页免| 成年午夜精品久久精品| 欧美日韩国产成人在线观看| 亚洲黄色视频在线观看一区| 成年A级毛片| 日韩精品成人网页视频在线| 99r在线精品视频在线播放| 麻豆精品在线播放| 欧美日韩国产在线观看一区二区三区| 国产福利在线免费| 伊人久久婷婷五月综合97色| 日韩天堂视频| 狠狠色香婷婷久久亚洲精品| 久久精品国产精品国产一区| 久久99精品久久久久久不卡| 久无码久无码av无码| 亚洲一级色| 2020久久国产综合精品swag| 久久国产精品国产自线拍| 国产三级国产精品国产普男人 | 国产精品美人久久久久久AV| 五月婷婷丁香色| 一级毛片在线播放免费| 欧美精品v欧洲精品| 亚洲人妖在线| 2020精品极品国产色在线观看 | 999精品视频在线| 乱人伦视频中文字幕在线| 国产99热| 不卡网亚洲无码| 欧美色香蕉| 国产SUV精品一区二区| www.99在线观看| 国产福利免费视频|