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

衛(wèi)星星座分布概率解析算法及精度分析

2023-11-18 05:24:02劉慧梁孫茜楚堯彭菲江帆呂紅劍蔡亞星王平
中國空間科學(xué)技術(shù) 2023年5期

劉慧梁,孫茜,楚堯,彭菲,江帆,呂紅劍,蔡亞星,王平

1.中國空間技術(shù)研究院 通信與導(dǎo)航衛(wèi)星總體部,北京 100094 2.國家航天局 衛(wèi)星通信系統(tǒng)創(chuàng)新中心,北京 100094 3.啟元實(shí)驗(yàn)室,北京 100095

1 引言

近50年來,人類探索太空的歷史先后經(jīng)歷了兩次低軌衛(wèi)星星座的浪潮。20世紀(jì)90年代,以“銥星”“全球星”為代表的第一代低軌星座系統(tǒng)旨在提供全球覆蓋的通信服務(wù),其星座規(guī)模在數(shù)十顆衛(wèi)星量級(jí),提供窄帶通信服務(wù),有效補(bǔ)充了地面蜂窩網(wǎng)絡(luò)未覆蓋地區(qū)的通信服務(wù)空白,特別是為南北極地區(qū),提供了高效、便捷的通信服務(wù)[1]。但由于終端成本及業(yè)務(wù)運(yùn)營等問題,上述衛(wèi)星通信公司先后都經(jīng)歷了破產(chǎn)重組,低軌衛(wèi)星星座系統(tǒng)的第一次發(fā)展浪潮也隨之陷入低谷。2015年,在谷歌等互聯(lián)網(wǎng)巨頭的推動(dòng)和支持下,以O(shè)neWeb和Starlink為代表的低軌星座在短期內(nèi)迅速聚集人氣,掀起了低軌衛(wèi)星星座開發(fā)的第二波浪潮。與第一代低軌星座系統(tǒng)相比,第二代低軌衛(wèi)星星座規(guī)模有了大幅提升,多數(shù)是由數(shù)千顆甚至上萬顆衛(wèi)星組成的大型星座,主要提供寬帶互聯(lián)網(wǎng)接入服務(wù)[2-4]。

伴隨著國外大型低軌星座的快速部署,近地軌道將面臨多系統(tǒng)共存的局面[5],為了避免頻率軌道資源干擾沖突等諸多問題,通常在星座設(shè)計(jì)階段會(huì)對在軌及規(guī)劃的星座進(jìn)行分析,根據(jù)星座特征參數(shù)構(gòu)建分析模型,其中,星座系統(tǒng)衛(wèi)星位置分布概率的預(yù)測及其準(zhǔn)確度直接影響著衛(wèi)星互聯(lián)網(wǎng)星座系統(tǒng)的運(yùn)行策略規(guī)劃[6-7]、通信協(xié)議設(shè)計(jì)[8-9]以及頻率干擾分析[10-11]等系統(tǒng)級(jí)設(shè)計(jì)驗(yàn)證工作。高效、準(zhǔn)確地實(shí)現(xiàn)多系統(tǒng)間衛(wèi)星位置分布概率提取,是高質(zhì)量完成系統(tǒng)設(shè)計(jì)的基礎(chǔ),具體分析方法可以分為基于軌道外推的時(shí)域統(tǒng)計(jì)法以及基于位置概率的解析算法。

基于軌道外推法的時(shí)域統(tǒng)計(jì)法通過考慮衛(wèi)星的各種攝動(dòng)模型,建立相應(yīng)的軌道方程,精度高,但計(jì)算復(fù)雜度大,對于目前動(dòng)輒萬顆衛(wèi)星規(guī)模的大型低軌星座來說,時(shí)域統(tǒng)計(jì)法計(jì)算耗時(shí)長、對服務(wù)器算力要求高,無法進(jìn)一步支持對星座系統(tǒng)覆蓋率、通信協(xié)議、星座頻率干擾等開展分析[12-14]。

為彌補(bǔ)上述時(shí)域統(tǒng)計(jì)方法的不足,提高星座分布概率預(yù)測的效率,國際電聯(lián)先后研究并頒布了S.1529建議書“非靜止軌道(non-geostationary orbit,NGSO)衛(wèi)星固定業(yè)務(wù)系統(tǒng)和其他衛(wèi)星系統(tǒng)間干擾統(tǒng)計(jì)測定的分析方法”[15],S.1257建議書“從地球表面看NGSO衛(wèi)星的可見度及干擾數(shù)據(jù)的計(jì)算分析方法”[16],給出了一種基于位置概率的計(jì)算方法。通過對典型場景的驗(yàn)證,基于位置概率的計(jì)算方法與時(shí)域統(tǒng)計(jì)方法的結(jié)果基本一致,計(jì)算時(shí)間大幅降低[13,17]。然而,國際電聯(lián)相關(guān)研究在驗(yàn)證基于位置概率的計(jì)算結(jié)果時(shí),所用時(shí)域統(tǒng)計(jì)仿真結(jié)果包含大量假設(shè)前提,例如不考慮地球自轉(zhuǎn)的影響(即假設(shè)地球自轉(zhuǎn)角速度為0),設(shè)定升交點(diǎn)漂移為固定值(0.06°每公轉(zhuǎn)周期)等,與現(xiàn)實(shí)場景之間存在較大差距。此外,兩份建議書只對特定場景的概率計(jì)算結(jié)果與時(shí)域統(tǒng)計(jì)法進(jìn)行了數(shù)值對比,并未對計(jì)算精度及誤差產(chǎn)生的原因做系統(tǒng)分析,所提方法對由數(shù)千顆衛(wèi)星構(gòu)成的大型星座的適用性也沒有驗(yàn)證。

本文參考國際電聯(lián)前期所開展的基于衛(wèi)星星座位置概率的研究基礎(chǔ),首次提出并詳細(xì)推導(dǎo)了具有普適性的衛(wèi)星星座概率分布的解析表達(dá)式,引入了地面站可視空域小區(qū)內(nèi)衛(wèi)星期望值變量,系統(tǒng)驗(yàn)證了對于不同構(gòu)型大型星座的適用性,通過與時(shí)域統(tǒng)計(jì)法結(jié)果對比,詳細(xì)分析了地面站在不同緯度下面對單一構(gòu)型星座場景及混合構(gòu)型星座場景的計(jì)算精度,為精確解析大型衛(wèi)星星座分布概率提供了理論基礎(chǔ)及高效計(jì)算方法。

2 解析法公式推導(dǎo)

在衛(wèi)星系統(tǒng)間頻率干擾分析中,利用衛(wèi)星分布概率確定星座空間位置關(guān)系,模擬干擾概率分布曲線,是完成星座系統(tǒng)性能優(yōu)化及干擾規(guī)避策略設(shè)計(jì)一種重要方法[13]。本節(jié)將詳細(xì)推導(dǎo)解析法求解不同構(gòu)型星座衛(wèi)星在可視空域小區(qū)中的分布特征,首先求解衛(wèi)星星座概率密度函數(shù),之后在此基礎(chǔ)上推導(dǎo)地面站可視空域小區(qū)內(nèi)衛(wèi)星期望值。

圖1為衛(wèi)星地面站可視空域劃分示意圖。將任意地面站可視空域劃為更小的小區(qū)單元,以正六邊形為網(wǎng)格劃分,小區(qū)單元緊密排列,各區(qū)域中心呈三角形排布,方位角0°表示正北方向,90°表示正東方向,180°表示正南方向,極軸長度表示仰角范圍,坐標(biāo)中心點(diǎn)表示地面站正上方,即仰角為90°。

圖1 衛(wèi)星地面站可視空域劃分示意Fig.1 Schematic of the Field of View(FOV) division of the satellite earth terminal

2.1 衛(wèi)星星座概率密度函數(shù)

假設(shè)某衛(wèi)星運(yùn)行在軌道半徑為R、傾角為δ的圓形軌道上,若衛(wèi)星軌道為非回歸軌道,則衛(wèi)星軌跡將逐步遍歷半徑為R的截?cái)嗲蛎?形成軌道殼,如圖2所示。

圖2 非靜止軌道衛(wèi)星軌道殼示意Fig.2 Schematic of the orbit shell for NGSO (non-geostationary orbit)satellite

該軌道殼上的衛(wèi)星位置可用星下點(diǎn)經(jīng)度ψ和緯度φ表示,由于圍繞地球極軸的自轉(zhuǎn)具有對稱性,所以衛(wèi)星軌道殼的經(jīng)度將均勻分布,但是,衛(wèi)星所處緯度的概率,與其運(yùn)行軌道傾角相關(guān)。運(yùn)行于圓軌道的衛(wèi)星位置矢量如下:

式中:τ為服從[0,1]均勻分布的隨機(jī)變量;z為衛(wèi)星位置矢量在地心指向北極方向上的分量。

(1)

式(1)存在兩個(gè)解,分別對應(yīng)運(yùn)行中衛(wèi)星星下點(diǎn)緯度增大和緯度減小的兩種情況,因此衛(wèi)星位于緯度區(qū)間[φ1,φ2]內(nèi)的概率為:

則衛(wèi)星在緯度方向的概率密度為:

由于衛(wèi)星在軌道殼經(jīng)度方向均勻分布,所以,概率密度函數(shù)可表示為:

(2)

圖3展示了傾角為45°的圓形軌道面上衛(wèi)星的概率密度函數(shù),其在南北緯大于45°區(qū)域概率密度為0;在小于或等于45°的區(qū)域,赤道上空衛(wèi)星概率密度最低,隨緯度增加,概率密度單調(diào)遞增。

圖3 傾角為45°時(shí)衛(wèi)星概率密度函數(shù)Fig.3 Probability density function of a satellite placed in a circular orbit plane with the inclination of 45°

2.2 地面站可視空域小區(qū)內(nèi)衛(wèi)星期望值

地面站T與其可視空域小區(qū)Ci的幾何關(guān)系如圖4所示。衛(wèi)星出現(xiàn)在小區(qū)Ci中的概率為

圖4 地面站可視空域小區(qū)幾何關(guān)系示意Fig.4 Schematic of geometric relationship for NGSO satellite in the FOV of satellite earth terminal

p=pC(ψ,φ)dS

(3)

即概率密度函數(shù)在小區(qū)Ci上的面積分。假設(shè)概率密度函數(shù)在被積分面積上均勻,則有

pi=pC(ψi,φi)AC

(4)

式中:pC(ψi,φi)為小區(qū)單元中心點(diǎn)對應(yīng)的衛(wèi)星概率密度;AC為地面站可見空域小區(qū)單元的球面度(單位sr)。若可視空域小區(qū)為圓形,則面積AC的表達(dá)式為:

式中:Δθε和Δθβ分別為空域小區(qū)Ci的俯仰角差和方位角差所對應(yīng)的地心角差,可由地面站與可見空域小區(qū)的位置幾何關(guān)系求出,表達(dá)式如下:

θε_(tái)min=arccos[kcos(ε-rcell)]-(ε-rcell)

θε_(tái)max=arccos [kcos(ε+rcell)]-(ε+rcell)

Δθε=θε_(tái)max-θε_(tái)min

式中:r與R分別為地球半徑和衛(wèi)星圓軌道半徑,ε和θε分別為空域小區(qū)中心點(diǎn)對應(yīng)的地面站仰角與地心角,θε_(tái)min與θε_(tái)max分別為空域小區(qū)上邊緣與下邊緣對應(yīng)的地心角,rcell為可視空域圓形小區(qū)半徑。

當(dāng)衛(wèi)星總數(shù)為N時(shí),對某一確定地面站,其可視空域小區(qū)Ci內(nèi),衛(wèi)星期望值為:

E=Npi=NpC(ψi,φi)AC

(5)

當(dāng)衛(wèi)星星座存在k種不同構(gòu)型時(shí),對某一確定地面站,其可視空域小區(qū)Ci內(nèi),衛(wèi)星期望值為:

(6)

3 不同構(gòu)型星座計(jì)算結(jié)果及精度分析

本節(jié)將利用以上推導(dǎo)的解析法公式,求解不同構(gòu)型星座對于地面站可視空域小區(qū)衛(wèi)星出現(xiàn)期望值的計(jì)算結(jié)果,并對計(jì)算精度進(jìn)行分析。本文選取兩個(gè)典型的Walker星座構(gòu)型作為研究對象,分別研究單一構(gòu)型及混合構(gòu)型條件下的分析結(jié)果,具體構(gòu)型軌道參數(shù)見表1。其中,構(gòu)型A為δ型Walker星座,包括3200顆衛(wèi)星,軌道高度1150km,傾角為60°;構(gòu)型B為星型Walker星座,包括2880顆衛(wèi)星,軌道高度1200km,傾角為88°。

表1 非靜止軌道通信星座軌道參數(shù)Table 1 Orbit parameters of NGSO satellite systems

在驗(yàn)證計(jì)算精度時(shí),采用基于軌道外推的時(shí)域統(tǒng)計(jì)結(jié)果ES與基于位置概率的解析法計(jì)算結(jié)果EA的比值作為評估指標(biāo),其中基于軌道外推的時(shí)域統(tǒng)計(jì)結(jié)果采用圓軌道,仿真步長為10s,總仿真時(shí)長為10000h(約417d),即以年為星座運(yùn)行時(shí)長量級(jí)評估解析算法計(jì)算精度。在地面站空域劃分中,本文選取可視空域小區(qū)中心點(diǎn)仰角大于20°區(qū)域范圍進(jìn)行統(tǒng)計(jì)及計(jì)算,可視空域小區(qū)半徑rcell最大選取10°,最小選取2°。圖5(a)展示了可視空域小區(qū)半徑為10°時(shí),仰角大于20°區(qū)域范圍被劃分61個(gè)小區(qū)單元,圖5(b)展示了可視空域小區(qū)半徑為4°時(shí),相同區(qū)域范圍被劃分367個(gè)小區(qū)單元。小區(qū)單元緊密排列,各區(qū)域中心呈三角形排布,極坐標(biāo)中,極軸長度表示仰角范圍,坐標(biāo)中心點(diǎn)表示地面站正上方,即仰角為90°,方位角0°表示正北方向,90°表示正東方向,180°表示正南方向。

圖5 地面站可視空域小區(qū)單元?jiǎng)澐諪ig.5 Schematic of the FOV division of the satellite earth terminal

3.1 單一構(gòu)型星座計(jì)算結(jié)果及精度分析

首先,對單一構(gòu)型星座衛(wèi)星位置概率計(jì)算結(jié)果進(jìn)行分析。構(gòu)型A星座在不同緯度地面站可視空域衛(wèi)星期望值由式(5)計(jì)算,具體期望值分布如圖6所示。結(jié)果中圓圈位置表示子空域中心,灰度圖例表示衛(wèi)星出現(xiàn)在小區(qū)單元的期望值,方位角0°表示正北方向,極軸長度表示仰角范圍。當(dāng)?shù)孛嬲疚挥?°(N)時(shí),正上方空域(仰角等于90°)小區(qū)單元內(nèi)衛(wèi)星期望值最小,隨著仰角降低,四周小區(qū)單元衛(wèi)星期望值逐漸增加,由于構(gòu)型A星座包含衛(wèi)星數(shù)量高達(dá)3200顆,所以在視場邊緣小區(qū)半徑為10°區(qū)域(正北、正南方向)衛(wèi)星期望值可達(dá)3.2。當(dāng)?shù)孛嬲疚挥?0°(N)時(shí),隨著地面站緯度增加,圖中北方的衛(wèi)星出現(xiàn)概率會(huì)略高于南方。當(dāng)?shù)孛嬲疚挥?0°(N)時(shí),衛(wèi)星集中于可視范圍內(nèi)南方區(qū)域,北方大部分區(qū)域衛(wèi)星期望值為0,對應(yīng)式(2)中衛(wèi)星星下點(diǎn)緯度大于衛(wèi)星軌道傾角的區(qū)域概率密度為0。

圖6 構(gòu)型A星座地面站可視空域小區(qū)衛(wèi)星期望值Fig.6 Expectation of satellite in the FOV of satellite earth terminal for the constellation with configuration A

圖7展示了采用基于軌道外推的時(shí)域統(tǒng)計(jì)小區(qū)單元內(nèi)衛(wèi)星期望值ES與基于位置概率的解析法計(jì)算衛(wèi)星期望值EA的比值,可以看出,小區(qū)單元半徑越小,則期望值比值越接近于1,即基于軌道外推的時(shí)域統(tǒng)計(jì)結(jié)果與基于衛(wèi)星位置概率的解析算法越接近。從理論上講,即小區(qū)單元面積越小,基于式(4)的近似算法與式(3)的概率密度面積分結(jié)果越吻合。地面站緯度小于或等于45°時(shí),基于位置概率的解析算法與基于軌道外推的時(shí)域統(tǒng)計(jì)結(jié)果相比,誤差最小為0.1%(小區(qū)半徑為2°時(shí)),誤差最大為3.1%(小區(qū)半徑為10°時(shí))。當(dāng)?shù)孛嬲揪暥鹊扔谛l(wèi)星軌道傾角,即60°時(shí),由于衛(wèi)星軌跡覆蓋邊緣存在未占滿小區(qū)單元,總體誤差平均值有所增加,為4.1%至6.7%。

圖7 構(gòu)型A星座地面可視空域衛(wèi)星期望值基于軌道外推法計(jì)算結(jié)果ES與解析法計(jì)算結(jié)果EA比值Fig.7 Ratio of expectation by simulation and expectation by analytical method for constellation with configuration A

構(gòu)型B星座在不同緯度地面站可視空域衛(wèi)星期望值同樣由式(5)計(jì)算,具體期望值分布如圖8所示。當(dāng)?shù)孛嬲疚挥?°(N)時(shí),正上方空域衛(wèi)星出現(xiàn)概率最低,隨著仰角降低,四周小區(qū)單元內(nèi)衛(wèi)星出現(xiàn)概率逐漸增加,構(gòu)型B星座包含衛(wèi)星數(shù)量為2880顆,略小于構(gòu)型A星座,視場邊緣小區(qū)半徑為10°區(qū)域(正北、正南方向)衛(wèi)星期望值最大為2.6。當(dāng)?shù)孛嬲疚挥?0°(N)時(shí)可明顯看出,隨著地面站緯度增加,圖中北方的衛(wèi)星期望值會(huì)逐漸高于南方。由于構(gòu)型B星座傾角為88°,空間分布特征和之前所述構(gòu)型A星座表現(xiàn)出一定差異,其在高緯度地區(qū)仍有較好的可見性,當(dāng)?shù)孛嬲疚挥?0°(N)時(shí),正北方向可視空域邊緣小區(qū)半徑為10°區(qū)域衛(wèi)星期望值可達(dá)到12.2。

圖8 構(gòu)型B星座地面站可視空域小區(qū)衛(wèi)星期望值Fig.8 Expectation of satellite in the FOV of satellite earth terminal for the constellation with configuration B

圖9展示了構(gòu)型B星座采用基于軌道外推的時(shí)域統(tǒng)計(jì)小區(qū)單元內(nèi)衛(wèi)星期望值ES與基于位置概率的解析法計(jì)算衛(wèi)星期望值EA的比值。當(dāng)?shù)孛嬲揪暥葹?°(N)和30°(N),誤差最小為0.1%(小區(qū)半徑為2°時(shí)),誤差最大為2.9%(小區(qū)半徑為10°時(shí))。當(dāng)?shù)孛嬲揪暥葹?0°(N),誤差最小為0.05%(小區(qū)半徑為2°時(shí)),誤差最大為2.6%(小區(qū)半徑為10°時(shí))。

圖9 構(gòu)型B星座地面可視空域衛(wèi)星期望值基于軌道外推法計(jì)算結(jié)果ES與解析法計(jì)算結(jié)果EA比值Fig.9 Ratio of expectation by simulation and expectation by analytical method for constellation with configuration B

3.2 混合構(gòu)型星座解析精度分析

在對單一構(gòu)型星座衛(wèi)星位置概率計(jì)算結(jié)果

進(jìn)行分析的基礎(chǔ)上,本小節(jié)將針對包含構(gòu)型A及構(gòu)型B的混合構(gòu)型星座,展示由式(6)計(jì)算的不同緯度地面站可視空域衛(wèi)星期望值,并對計(jì)算精度進(jìn)行分析。

混合構(gòu)型星座期望值分布如圖10所示。由于該混合構(gòu)型星座為構(gòu)型A星座與構(gòu)型B星座疊加,所以圖10地面站可視空域小區(qū)衛(wèi)星期望值為圖6與圖8結(jié)果疊加。當(dāng)?shù)孛嬲疚挥?°(N)時(shí),正上方空域(仰角等于90°)小區(qū)單元內(nèi)衛(wèi)星期望值最小,隨著仰角降低,四周小區(qū)單元衛(wèi)星期望值逐漸增加,在視場邊緣小區(qū)半徑為10°區(qū)域(正北、正南方向)衛(wèi)星期望值為5.8。當(dāng)?shù)孛嬲疚挥?0°(N)時(shí),圖中北方的衛(wèi)星出現(xiàn)概率會(huì)略高于南方。當(dāng)?shù)孛嬲疚挥?0°(N)時(shí),衛(wèi)星期望值分布包括構(gòu)型A特征,即衛(wèi)星分布集中于可視范圍內(nèi)南方區(qū)域,又包括構(gòu)型B特征,正北方向可視空域邊緣小區(qū)半徑為10°區(qū)域衛(wèi)星期望值可達(dá)到12.2。

圖10 混合構(gòu)型星座地面站可視空域小區(qū)衛(wèi)星期望值Fig.10 Expectation of satellite in the FOV of satellite earth terminal for the constellation with hybrid configuration

圖11展示了混合構(gòu)型星座采用基于軌道外推的時(shí)域統(tǒng)計(jì)小區(qū)單元內(nèi)衛(wèi)星期望值ES與基于位置概率的解析法計(jì)算衛(wèi)星期望值EA的比值。當(dāng)?shù)孛嬲疚挥?°(N)和30°(N)時(shí),基于位置概率的解析算法誤差最小為0.1%(小區(qū)半徑為2°時(shí)),誤差最大為2.9%(小區(qū)半徑為10°時(shí))。當(dāng)?shù)孛嬲揪暥鹊扔跇?gòu)型A衛(wèi)星軌道傾角,即60°時(shí),總體誤差平均值增加為3.5%至5.8%。

圖11 混合構(gòu)型星座地面可視空域衛(wèi)星期望值基于軌道外推法計(jì)算結(jié)果ES與解析法計(jì)算結(jié)果EA比值Fig.11 Ratio of expectation by simulation and expectation by analytical method for constellation with hybrid configuration

4 結(jié)論

為滿足日益增長的衛(wèi)星星座空間分布概率快速分析需求,本文首次提出并詳細(xì)推導(dǎo)了具有普適性的星座概率分布解析表達(dá)式,引入了地面站可視空域小區(qū)內(nèi)衛(wèi)星期望值變量,并系統(tǒng)驗(yàn)證了所提解析算法對于不同構(gòu)型大型星座的適用性。與傳統(tǒng)基于軌道外推的時(shí)域統(tǒng)計(jì)法結(jié)果相比,期望值解析算法的結(jié)果誤差大小與可視空域小區(qū)單元半徑相關(guān),小區(qū)單元半徑越小,解析算法結(jié)果誤差越小:當(dāng)可視空域小區(qū)單元半徑等于10°時(shí),最大誤差為3.1%,當(dāng)小區(qū)半徑減少至2°時(shí),最大誤差小于0.1%。除此之外,本文還詳細(xì)分析了地面站部署在不同緯度時(shí)對單一構(gòu)型星座場景及混合構(gòu)型星座場景的計(jì)算精度的影響。當(dāng)?shù)孛嬲揪暥葟?°不斷增加至星座軌道傾角角度時(shí),衛(wèi)星軌跡覆蓋邊緣會(huì)逐漸出現(xiàn)未占滿小區(qū)單元、從而導(dǎo)致計(jì)算結(jié)果誤差增加的情況。利用本文所提解析算法可以快速計(jì)算地面站可視空域內(nèi)特定區(qū)域衛(wèi)星出現(xiàn)期望值,為快速評估衛(wèi)星星座分布概率、星座系統(tǒng)覆蓋性分析及系統(tǒng)間頻率干擾分析提供了理論基礎(chǔ)及高效計(jì)算方法。

主站蜘蛛池模板: 亚洲国产日韩在线观看| 欧美区一区| 欧美亚洲欧美| 免费高清a毛片| 99久久精彩视频| 国产精品浪潮Av| 亚洲三级网站| 国产精品久久国产精麻豆99网站| 国内精品视频| 5555国产在线观看| 亚洲色无码专线精品观看| 成人在线视频一区| 国产成人a在线观看视频| 国产在线视频二区| 一本大道视频精品人妻 | 风韵丰满熟妇啪啪区老熟熟女| 精品自拍视频在线观看| 国产精品主播| 黄色网站不卡无码| 亚洲资源站av无码网址| 国产乱论视频| 欧美一级专区免费大片| 91视频日本| 欧美另类图片视频无弹跳第一页 | 毛片久久网站小视频| 亚洲天堂.com| 精品国产黑色丝袜高跟鞋| 国产成人91精品免费网址在线 | 久久久久久久97| 欧美全免费aaaaaa特黄在线| 男人天堂亚洲天堂| 国产午夜无码专区喷水| 亚洲视频四区| 园内精品自拍视频在线播放| 欧美三級片黃色三級片黃色1| 91久久精品日日躁夜夜躁欧美| 久久久久久高潮白浆| 国产欧美日韩资源在线观看| 久久国产精品77777| 免费播放毛片| 亚洲经典在线中文字幕| 中文字幕乱妇无码AV在线| 久久99热这里只有精品免费看| 国产97视频在线| 亚洲成人黄色网址| 欧美国产视频| 国产精品成人第一区| 亚洲大尺码专区影院| 国产主播一区二区三区| 在线免费观看AV| 亚洲三级影院| 午夜不卡福利| 91在线精品麻豆欧美在线| 五月婷婷综合网| 亚洲最猛黑人xxxx黑人猛交| 午夜三级在线| 国产激情无码一区二区免费| 9久久伊人精品综合| 欧洲一区二区三区无码| 久爱午夜精品免费视频| 色哟哟国产精品| 久久精品91麻豆| 538国产在线| 国产成人精品一区二区| 久久香蕉国产线看观看精品蕉| 久久99久久无码毛片一区二区 | 怡红院美国分院一区二区| 91毛片网| 欧美亚洲激情| 国产精品色婷婷在线观看| 一级看片免费视频| 国产乱人免费视频| 国产精品视频第一专区| 香蕉国产精品视频| 草草影院国产第一页| 亚洲国产亚综合在线区| 日韩不卡高清视频| a级毛片免费看| 99re这里只有国产中文精品国产精品| 色AV色 综合网站| 精品少妇人妻一区二区| 欧美、日韩、国产综合一区|