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

基于LSTM神經網絡的短期軌道預報

2022-03-11 01:50:50張心宇宋佳凝
系統工程與電子技術 2022年3期
關鍵詞:模型

張心宇, 劉 源,*, 宋佳凝

(1. 中山大學天琴中心, 廣東 珠海 519082; 2. 中山大學物理與天文學院, 廣東 珠海 519082)

0 引 言

近年來隨著國內外航天技術的發展,空間任務、空間系統均日趨復雜化,這對衛星在軌的智能性、自主性提出了越來越高的要求。為了提升衛星的在軌自主性,使其具備較強的在軌自主任務規劃與任務協同能力,則衛星至少需要具備中短期的在軌自主軌道預報能力。而現有軌道預報方法主要存在預報精度低與計算量大的缺陷[1-4]。

目前傳統軌道預測方法主要是基于物理建模,一般有兩種思路[1]。一種是數值法[2],從初始狀態開始,利用狀態變化的導數逐步積分,從而得到后來任意時刻的狀態。利用數值法進行高精度軌道預報不僅需要對衛星幾何信息、地心引力,大氣阻力、太陽輻射壓力等方面有精細的模型,更會消耗巨大的計算成本和時間成本。另一種是解析法[3],根據影響軌道運動的主要攝動力得出衛星運動的解析表達式,通過求解析解對衛星軌道進行預報。解析法的優勢在于計算速度極快,運行時只占用少量的計算資源。但是由于各種攝動力的精確模型十分復雜,準確求解出解析解極其困難,只能對攝動模型進行一定程度上的簡化,得到精度較低的近似解,更多被應用在數量巨大、預報精度要求低的空間碎片預報上。因此,無論是數值法與解析法均難以滿足衛星平臺自主高精度軌道預報的需求[4]。神經網絡作為一種機器學習算法,能夠在不進行具體建模的基礎上對復雜規律進行學習,掌握其輸入輸出函數的內在聯系。同時,神經網絡在進行訓練時計算量較大,但是完成訓練后利用其進行計算時計算量則較小,注入衛星后能夠比現在主流的四階Runge-Kutta數值法與解析法方法更適用于星上的計算環境[4]。

Peng等利用人工神經網絡學習傳統方法預測衛星軌道的歷史誤差,顯著提升了衛星預測精度[5]。除了神經網絡,他們還利用支持向量機(support vector machine, SVM)[6-8]以及高斯過程(Gaussian process, GP)[9]兩種機器學習算法對傳統方法的歷史預測誤差進行學習,并將3種機器學習算法的結果進行了詳細的對比[10]與缺陷分析[11]。類似的還有Li等人提出的一種基于機器學習的單站稀疏跟蹤數據改進低軌空間碎片軌道預測的方法[12]。文獻[13]提出基于補償波形調整的導航衛星軌道預報方法,采用補償波形調整的神經網絡優化模型后,預報弧長為8 d、15 d、30 d時,改進率分別提高了2.3%、6.7%、10%。曹磊[14]以全球定位系統(global positioning system, GPS)導航衛星為對象,分別采用反向傳播(back propagation, BP)網絡和深度神經網絡對預報誤差進行建模并進行了預報,結果表明基于深度神經網混合模型的方法對于GPS導航衛星的預報精度有較大程度的改進。Yang在沒有物理建模的情況下,直接利用長短期記憶(long short-term memory, LSTM)學習衛星軌道變化規律,并用螢火蟲算法加以改進,一天預報精度達到了百米量級[15]。文獻[16]提出將優化的LSTM網絡直接應用在衛星軌道數據上,將預報20 d的誤差由之前最大值的300 km降低到5 km以下。文獻[15-16]利用了LSTM神經網絡,但都沒有利用到任何動力學模型,直接利用預報軌道作為訓練樣本,而且除了軌道坐標數據沒有引入任何特征,因此導致神經網絡輸出量的動態范圍過大,預報精度仍有進一步提升的空間。同時,之前提出的基于神經網絡的軌道預報方法都需要長達幾百天訓練的數據,限制了模型應用的靈活性。針對上述不足,本文提出一種新的軌道預報算法,在SGP4(simplified general perturbations)解析法軌道預報模型的基礎上,利用LSTM神經網絡以位置、速度、加速度作為特征,對SGP4軌道預報歷史誤差進行學習,通過預測并修正未來SGP4軌道預報誤差,改進動力學模型的短期軌道預報精度。

1 誤差數據分析

1.1 SGP4預報誤差

北美防空司令部(North American Aerospace Defense Command, NORAD)提供的TLE(two-line mean element)是目前最為完整的地球軌道空間目標的軌道數據,SGP4/SDP4(simplified deep-space perturbation)模型[17]是NORAD提供的配合 TLE 使用的解析軌道預報模型,能夠快速計算出任何時刻衛星的狀態。由于模型計算效率較高,預報精度尚可,因此被廣泛應用[18]。SGP4/SDP4算法對于近地目標預報3 d位置誤差小于40 km;對半同步軌道預報30 d和同步軌道預報15 d,預報誤差小于40 km;橢圓軌道目標預報1 d位置誤差小于20 km[19]。

TLE+SGP4預報誤差主要由兩部分組成[20]。第一部分是TLE軌道根數的誤差,是由監視網的測量誤差與定軌中擬合算法誤差共同引起的[21]。NORAD并沒有提供更詳細關于TLE誤差的說明,在沒有更高精度測量數據的情況下,這部分誤差是無法消除的。第二部分是SGP4模型的誤差,主要由采用不精確的大氣阻力模型與低階的地球引力場模型導致[22]。SGP4模型為了確保較快的運算速度,采用了簡化的攝動力模型。其大氣模型只考慮了大氣阻尼的長期影響,采用大氣密度隨高度變化的指數大氣密度模型。模型中氣動阻力項B*為歸一化的大氣阻力系數[23]:

(1)

式中:

(2)

ρo為大氣密度參考值;B為彈道系數;Cd為無量綱的氣動阻力系數;A/m為衛星面質比。B*是利用歷史數據將其與其他軌道根數一同進行解算,獲得平均化的計算初值,在SGP4外推過程一直為不變項,這必然會引起預報的較大誤差[24]。

同時SGP4使用的地球引力模型只考慮了J2、J3、J4帶諧項的影響,并未考慮J22田諧項的影響。對于典型的102量級(歸一化量綱)的預報,由忽略田諧項產生的誤差可以達到百米級甚至千米級[25]。本文希望通過學習歷史誤差數據的變化規律,提前預測TLE+SGP4的預報誤差,從而提高短期軌道預報精度。

1.2 誤差時間序列數據

時間序列是一組按照時間先后順序排列且內部關聯的數據。通過對歷史時間序列的分析,能在一定統計意義上對未來進行短期預測。地心慣性坐標系(earth center inertial coordinates, ECI)表示的衛星X、Y、Z三軸的位置誤差數據,也可以被認為是3組連續的時間序列。

表1 Ajisai衛星相關參數

以Ajisai衛星為例,其相關參數如表1所示,該衛星的TLE數據每天會發布1~3個,使用SGP4模型對TLE數據進行預報,采用前一天最后一次發布的TLE根數進行外推可以通過SGP4模型對后一天做出精度最高的預報,最大預報誤差的量級在600 m,所以本文主要討論對外推一天的預報誤差進行修正。作為標準的高精度的軌道數據選擇國際激光測距服務(international laser ranging service,ILRS)提供的CPF(consolidated prediction format)星歷文件[26],可以從CDDIS(crustal dynamics data information system)或者EUROLAS數據中心得到。對于Ajisai衛星,每天更新的CPF星歷文件提供5 d的軌道數據,精度可以達到2 m以內[27],滿足本文訓練神經網絡需要的精度。一個特定的CPF文件包括未來3 d與過去兩天共5 d的精確軌道預測,步長為4 min,并以國際地面參考系表示。由于真實情況下未來時刻的精密數據是未知的,所以選擇連續多個CPF文件的第二天數據連接起來作為實驗的精確星歷。由連續7 d的CPF星歷與SGP4預報得到的X軸、Y軸、Z軸3組原始誤差時間序列數據,如圖1所示。本文將對圖中數據之后一天的SGP4+TLE預報誤差序列進行修正。

2 基于LSTM的軌道預報誤差修正

2.1 LSTM網絡模型

在循環神經網絡的基礎上,提出的LSTM網絡[28-29]采用巧妙的門設計,避開了梯度爆炸和長期依賴問題[30]。由于每一次循環都用到了前一次循環的信息,每一個輸出狀態都受到了之前狀態的影響,所以,LSTM網絡能夠更好地記住長期的規律,并廣泛應用于時間序列預測的問題,也適用于處理衛星軌道數據的誤差關系[15]。

本文采用的LSTM網絡結構如圖2所示,在t時刻,一層網絡由上下兩條信息流構成。從Ct-1到Ct的信息流表示細胞狀態的傳遞,整條線路通過3個門控結構與下面的一條信息流線性交互。門控結構讓信息進行選擇式通過,流入上方Ct-1到Ct的信息流對細胞狀態進行刪除或者添加信息。門控結構中σ激活函數層與tanh激活函數層,分別可以將輸入轉化到(0,1)之間與(-1,1)之間,生成輸入數據的權重,從而對輸入數據進行篩選。每一層LSTM網絡中有3個門結構來控制細胞狀態。

(1) 遺忘門

遺忘門的表達式如下所示:

ft=σ(Wf·[ht-1,xt]+bf)

(3)

σ層通過上一時刻隱藏狀態ht-1與t時刻輸入xt得到一個0到1之間的值,作為上一層細胞狀態被遺忘的概率,視為Ct-1記憶的衰減系數,記作ft。

(2) 輸入門

it=σ[Wi(ht-1,xt)+bi]

(4)

(5)

(6)

式中:*表示卷積。

(3) 輸出門

ot=σ[Wo(ht-1,xt)+bo]

(7)

ht=ot*tanh(Ct)

(8)

輸出門通過式(7)和式(8)更新當前時刻隱藏層輸出的值即ht,同時也是t時刻的隱藏狀態。

式(3)~式(8)中Wζ與bξ,ζ∈{f,i,c,o},ξ∈{f,i,c,o}分別表示輸出的權重和偏置矩陣,也是在訓練中需要學習的參數。由于每一次循環都用到了前一次循環的信息Ct-1與ht-1,每一個輸出狀態都受到了之前狀態的影響,所以LSTM具有記憶長期歷史信息的能力。本文的網絡模型主要由兩個LSTM層以及一個線性全連接層組成。

2.2 過擬合修正

網絡訓練過程,可以看作是不斷調整參數以使損失函數最小化的過程。假設樣本量為N,通過網絡正向傳播得到輸出y*,而期望輸出為y,則可以給出均方差損失函數:

(9)

但是僅使用這樣的損失函數最終很可能導致過擬合問題,即對訓練集數據效果很好,而對實際的測試集數據效果很差。因此,本文采用了帶正則化項的均方差損失函數,即

(10)

式中:ω為網絡中各項參數;λ為正則化項的權重。

2.3 特征選擇與數據預處理

區別于之前論文方法僅基于位置誤差數據,本文的LSTM網絡模型引入位置誤差、速度預測值與加速度預測值數據,實現從多個特征對原始序列數據進行分析預測,提高信息的維度。其中目標的位置速度預測值與加速度預測值數據依靠前一天的TLE數據配合SGP4模型得到。

數據的預處理主要包括軌道殘差數據的完整性和歸一化步驟。時間序列數據的完整性,用于確保每一時刻監測值的有效性,是開展后續研究工作的前提。CPF星歷提供4 min間隔的高精度數據,如果需要更高密度的軌道數據,根據官方文件[26],通過十點基線拉格朗日插值器可以在CPF文件覆蓋的任意時刻獲得衛星的狀態,并且選取兩個點中間時刻插值精度最高。數據的歸一化處理用于消除位置、速度、加速度數據不同量綱之間的影響,本文使用歸一化函數如下:

(11)

式中:Xorigin、min、max分別表示原始數據、訓練數據中最小值與最大值;Xm為歸一化后的值,處理后的數據將被限制在[0,1]范圍內,保證網絡訓練過程能夠快速收斂。

2.4 數據集結構設計

根據衛星3個軸的歷史軌道誤差數據建立模型預測軌道預報誤差,首先需要設計合理的數據集。由于3個軸的數據集結構相同,此處以ECI坐標系下X軸從UTC(universal time coordinated) 2020-04-06的00:00:00 到UTC 2020-04-13的24:00:00連續8 d數據設計訓練數據集與測試數據集,Y、Z軸的數據集構造工作不再贅述。

訓練數據集的結構如圖3所示。固定時間長度的真實歷史誤差數據作為訓練數據依次輸入進LSTM網絡,然后前向傳播對輸入序列的后一個時刻的誤差做出一步預測,最后通過預測值與真實誤差進行反向傳播更新網絡的權值。一組輸入時間序列與下一步預測的結果構成一步訓練數據。圖3中每一步的訓練數據結構如圖4所示。

t0時刻LSTM網絡通過輸出步長為n,由預報誤差、歷史速度預報值與歷史加速度預報值構成的向量大小為3n的輸入時間序列,進行前向傳播對下一個時間點t0+1即4 min后的預報誤差進行一步預測,得到t0+1時刻的預報誤差預測值VarP(t0+1)。而t0+1時刻的真實預報誤差Var1(t0+1)作為參考值,通過計算出預測值VarP(t0+1)的交叉熵損失函數進行反向傳播,對LSTM網絡的權重矩陣與偏置進行更新。完成一步訓練后,輸入序列與預測序列會向后一個時間點移動,直到完成整個訓練集數據。

使用3個獨立的LSTM網絡對訓練數據分別進行300次訓練,圖5展示了訓練后的3個LSTM網絡在所有3個軸的訓練數據上的性能,反映了模型對解析法軌道預報誤差的學習能力。

原位雜交結果顯示,在結腸癌組織中,miR-454-3p的表達水平顯著高于正常結腸組織。我們對原位雜交結果進行統計分析,分別取了5個高倍視野,計算了miR-454-3p陽性細胞數占總細胞數的比例,結果顯示,結腸癌組織中的miR-454-3p表達水平顯著高于正常結腸組織中的表達水平組間差異顯著(P<0.05,圖1)。

通過圖5的殘余誤差可以看出,LSTM在X軸訓練數據上的最大殘余誤差已經下降到30 m左右;在Y軸訓練數據上的最大殘余誤差已經下降到60 m左右;在Z軸訓練數據上的最大殘余誤差已經下降到60 m左右。

2.5 預報誤差的預測流程

以上分析表明,LSTM網絡可以較好地學習到訓練集上SGP4預報誤差的變化規律。為了驗證LSTM網絡對未來軌道預報誤差的修正能力,本文著重討論SGP4+TLE外推一天的軌道預報誤差修正,將訓練后的LSTM網絡應用在測試數據集上。

測試數據集的結構與圖3訓練數據結構類似,都是單步預測,整個預測過程不進行反向傳播,僅執行前向傳播過程。在t0時刻預測t0+1時刻誤差的公式為

VarP(t0+1)=LSTM{{Var1(t0),Var2(t0),Var3(t0)},
{Var1(t0+1),Var2(t0+1),Var3(t0+1)},…,
{Var1(t0-n+1),Var2(t0-n+1),Var3(t0-n+1)},
{Var1(t0-n),Var2(t0-n),Var3(t0-n)}}

(12)

其中,LSTM函數表示神經網絡整個前向傳播得到預測輸出的過程。

但是,由于未來時刻的SGP4模型真實預報誤差未知,因此UTC 2020-04-13的00:00:00以后時刻的真實預報誤差Var1是不可知的,而速度預報值Var2與加速度預報值Var3可以通過SGP4模型得到。UTC 2020-04-13的00:00:00預測點之后的誤差數據,需要通過LSTM的預測數據填充作為輸入數據,詳細預測流程如圖6所示。t0時刻LSTM網絡通過輸出步長為n,向量大小為3n的輸入時間序列,對下一個時間點t0+1即4 min后的預報誤差進行一步預測。預測公式即為式(12)。通過網絡前向傳播得到t0+1時刻的預報誤差預測值VarP(t0+1)。完成一次預報后,輸入序列窗口會后移動一個步長。同時,輸入中的未知參數Var1(t0+1)會由上一次預測值VarP(t0+1)填補。

VarP(t0+2)=LSTM{{VarP(t0+1),Var2(t0+1),

Var3(t0+1)},…,{Var1(t0-n+2),Var2(t0-n+2),

Var3(t0-n+2)},{Var1(t0-n+1),Var2(t0-n+1),

Var3(t0-n+1)}}

(13)

3 仿真實驗與分析

3.1 仿真實驗結果

為了評估神經網絡對預報誤差的修正能力,本文使用殘差率Pml(ξ)來量化LSTM神經網絡對誤差的修正能力,數學表達式如下:

(14)

定義式(14)為修正后殘差絕對值之和與原始誤差絕對值之和的比值。其中n是測試數據的大小,下標ξ∈{x,y,z}表示誤差的不同分量。Pml(ξ)值越小,表示經過訓練的LSTM模型的性能越好。

實驗選擇UTC 2020-04-06的00:00:00 到UTC 2020-04-12的24:00:00每4 min一組數據,共7 d 2 520組數據作為訓練數據(訓練集上的殘差如圖5所示);選擇UTC2020-04-13的00:00:00到UTC 2020-04-13的24:00:00每4 min一組數據,共360組數據作為測試數據。

通過訓練好的3個獨立的LSTM神經網絡對衛星UTC 2020-04-13的00:00:00以后400 min、800 min和1 d的SGP4軌道預報誤差進行預測。在圖7中,展示了3個經過訓練的LSTM模型在測試數據上的對誤差預測性能,紅色線為SGP4預報的真實誤差,藍色線為LSTM模型預測的預報誤差。圖7分別為對X、Y、Z軸進行時長400 min、800 min和1 440 min(即一天)的預測結果。對于不同預測時長,相應最優的LSTM網絡的輸入步長與神經網絡的神經元個數也不同,這部分會在下一節進行詳細說明。從圖7可以看出,經過訓練的LSTM模型可以較好地對ECI坐標系下3個軸的SGP4預報誤差進行預測。

在圖8中分別為對X、Y、Z軸進行了400 min、800 min、1 440 min(即一天)預測后的殘余誤差。X軸最大誤差分別下降到了50 m、75 m與100 m左右;Y軸最大誤差分別下降到了50 m、75 m與100 m左右;Z軸最大誤差分別下降到了50 m、75 m與100 m左右。

表2展示了在獨立進行10次訓練與預測并對實驗結果取均值的情況下,模型對不同時長誤差(網絡參數取值不一定是最優值)修正后的殘差率Pml(ξ)。結果表明,采用LSTM模型的修正方法可以大幅提高基于動力學方法的軌道預報精度。對X軸進行400 min、800 min、1 440 min的預報Pml(x)分別為10.26%、11.96%、16.87%。對Y軸進行400 min、800 min、1 440 min的預報Pml(y)分別為9.52%、13.25%、17.66%。對Z軸進行400 min、800 min、1 440 min的預報Pml(z)分別為9.30%、12.36%、19.58%。隨著預測時長的增加,測試數據的性能開始下降,這是可以預期的,因為訓練后的LSTM模型會使用預測的誤差值對后續的輸入序列中未知的Var1進行填補,導致誤差的迭代,從而使網絡性能下降。同時,預測步長與輸入序列步長的增加導致數據集樣本數量減少,同樣會導致性能下降。擴大訓練數據的數量通常可以提升模型性能,但是考慮到TLE數據與精密軌道文件每天都會更新的情況,以及盡量減少計算成本與時間成本的需求。在取舍下,本文選擇7的數據大小,需要的訓練集大小遠遠小于文獻[5-8]需要的100~250 d的訓練數據集。

表2 LSTM模型對不同時長誤差的預測能力

3.2 預報時長與輸入序列步長對LSTM模型性能影響

為了研究預報時長與輸入序列步長對模型性能的影響,本文以X軸數據為例,通過UTC 2020-04-06的00:00:00 到UTC 2020-04-12的24:00:00共7天2 520組數據作為訓練數據對LSTM神經網絡進行訓練并對衛星UTC 2020-04-13的00:00:00以后連續7 d即10 080 min的SGP4軌道預報的X軸誤差進行修正。

同時,通過改變輸入序列的步長研究輸入序列長度對LSTM模型性能影響。輸入序列步長即圖3與圖6中輸入窗口的長度,決定每次輸入進LSTM網絡中的數據大小。本文使用相同的X軸訓練集數據,同時分別選擇不同輸入序列長度的輸入對相同結構的LSTM網絡進行訓練,最后對UTC 2020-04-13的00:00:00以后連續7天數據進行誤差修正。上述實驗軌道預報修正后的殘余誤差如圖9所示,其中輸入序列步長分別為270、360、450。

通過圖9這3組的預測殘差可以看出,隨著預測時間增長模型修正后的殘余誤差逐漸增大。這是由于輸入窗口已經完全離開真實歷史誤差區域,輸入序列全部由之前的預測值填補,預測的誤差會進一步迭代,導致模型預測能力隨著誤差的累積降低,殘差率Pml(ξ)增大。同時輸入序列步長的增加可以提升模型的預測能力,當長度為270、360、450時,其Pml(ξ)分別為43.13%、30.45%與28.79%。

在表3中具體展示了不同輸入序列步長對殘差率Pml(x)的影響結果,可以看出模型性能隨著步長增加提升,印證了上述結論。

表3 不同輸入序列步長對殘差率Pml(x)的影響

4 結 論

為了解決現有的在軌自主高精度短期軌道預報精度較低、計算量較大的問題,本文給出了利用LSTM神經網絡對SGP4軌道預報進行修正的模型。將7 d的歷史軌道誤差、速度與加速度數據作為訓練的樣本,通過訓練3個獨立的LSTM網絡,對后一天的軌道數據誤差進行預測。經過預測值修正后ECI坐標系下3個軸的預報精度基本達100 m量級,殘余誤差Pml(ξ)下降到小于20%。同時,給出了預測一周的殘余誤差隨時間的變化趨勢以及輸入序列步長對殘差率Pml(x)的影響。

本研究表明,通過對歷史誤差規律進行學習并用以預測未來,能夠在原始動力學模型的基礎上大幅提高衛星軌道的預報精度。其預報精度基本滿足衛星自主軌道預報、地面測控需求的精度,同時又有訓練數據集更小、無需數值積分,適用于星上計算環境等優點,因此在應用上有著一定的參考價值。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 久久国产毛片| 亚洲天堂网2014| 国产无码精品在线| 亚洲天堂日韩av电影| 欧美中文字幕在线视频| 欧洲欧美人成免费全部视频| 国产69囗曝护士吞精在线视频| 亚洲人精品亚洲人成在线| 亚洲人成电影在线播放| 亚洲国产精品日韩av专区| 亚洲欧美国产五月天综合| 国产真实乱子伦精品视手机观看 | 亚洲一区二区精品无码久久久| 国产精品成人观看视频国产| 亚洲swag精品自拍一区| 爽爽影院十八禁在线观看| 亚洲欧美在线看片AI| 天天躁夜夜躁狠狠躁图片| 亚洲国产成人久久精品软件| 无码AV动漫| 欧美黄网在线| 亚洲精品国产精品乱码不卞| 免费一极毛片| 国产精品粉嫩| 欧美亚洲欧美| 中文字幕 91| 四虎亚洲国产成人久久精品| 成人精品区| 国产免费黄| 国产va在线观看| 国产精品久久精品| 精品99在线观看| 久久精品人妻中文视频| 国产亚洲欧美日韩在线观看一区二区| 五月婷婷导航| 色屁屁一区二区三区视频国产| 又爽又黄又无遮挡网站| 国产美女丝袜高潮| 亚洲成A人V欧美综合| 免费中文字幕在在线不卡| 亚洲视频欧美不卡| 91亚洲精品第一| 99re在线免费视频| 日本精品αv中文字幕| 色天堂无毒不卡| 潮喷在线无码白浆| 亚洲综合九九| 亚国产欧美在线人成| 亚洲免费毛片| 人妻精品久久无码区| 国产凹凸视频在线观看| 国产精品自拍露脸视频| 青青操国产| 国产免费a级片| 国产99在线| 国产精品欧美日本韩免费一区二区三区不卡| 亚洲另类色| 无码AV高清毛片中国一级毛片| 欧美 亚洲 日韩 国产| 97国产在线观看| 人人爱天天做夜夜爽| 狠狠ⅴ日韩v欧美v天堂| 国禁国产you女视频网站| 亚洲aaa视频| 欧美日韩专区| 男人天堂亚洲天堂| 国产精品女同一区三区五区| 国产另类视频| 欧美a在线| 免费一级毛片在线观看| 91久久天天躁狠狠躁夜夜| 欧美综合成人| 丝袜亚洲综合| 日韩AV无码免费一二三区| 欧美在线精品怡红院| 亚洲va视频| 国产高清在线精品一区二区三区 | 国产一级毛片在线| 中文字幕人妻无码系列第三区| 亚洲av日韩av制服丝袜| 国产精品黄色片| 日韩专区欧美|