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

基于非平穩時序ARIMA模型的W頻段雨衰預測

2015-10-14 04:05:03
電子與信息學報 2015年10期
關鍵詞:模型

楊 峰 薛 斌 劉 劍

?

基于非平穩時序ARIMA模型的W頻段雨衰預測

楊 峰 薛 斌*劉 劍

(空軍工程大學信息與導航學院 西安 710077)

針對目前絕大多數雨衰預測模型僅驗證到55 GHz,而經過驗證的W頻段預測模型相對較少,且存在模型表述復雜度高、計算量大的問題,該文提出一種結構簡單、計算量小的實時預測方法。該方法基于ARIMA模型,利用非平穩雨衰時序中相鄰時序間的相關性建立預測模型,對初始序列進行平穩性檢驗,通過差分變換將非平穩序列轉化為平穩序列,并對平穩化后的時間序列進行參數估計及診斷檢驗,將傳統非線性預測轉化為線性預測。并先將該ARIMA(1,1,6)模型在不同極化方式、預測間隔和時序個數的條件下進行比較,然后分別與ITU-R, Silva Mello模型在垂直極化、預測間隔0.10 GHz,時序個數50的條件下進行比較,最后使用ARIMA(1,1,6)模型進行預測,并對照預測序列與仿真序列的吻合度。結果表明,ARIMA模型與ITU-R, Silva Mello模型所得結果預測誤差不超過,且衰減變化趨勢基本相同,預測序列與仿真序列間吻合度較高,說明該方法可用于W頻段雨衰預測,且預測精度高,模型表述簡單。

雨衰預測;W頻段;非平穩時間序列;ARIMA;差分

1 引言

衛星通信已成為全球通信的重要組成部分,近年來,隨著寬帶多媒體衛星通信的蓬勃發展,C, Ku, Ka等較低頻段頻帶資源日益緊張。為了滿足衛星系統對通信容量和傳輸速率日益增長的需求,各國加強了對更高頻段的研究,如Q/V頻段和W頻段,目前W頻段已經成為衛星通信研究前沿[1]。美國,意大利等國已率先對W頻段衛星通信予以重點研究[1,2]。國際電信聯盟(International Telecommu- nication Union, ITU)為W頻段衛星通信分配的頻段是:上行81~86 GHz,下行71~76 GHz。與其他較低頻段相比,該頻段可用帶寬大,傳輸速率高,天線尺寸小,波束窄,抗干擾能力強[3,4],受到了各國廣泛重視。但與此同時它也存在一些問題,如由大氣吸收、雨衰、快衰落、反射等引起的惡劣的信道傳播特性[5],非線性失真、相位噪聲等嚴重的硬件問題[6]。其中降雨是影響通信質量最嚴重的大氣因素之一(數據體現,與其他頻段相比較),所以,準確預測降雨衰減對W頻段衛星通信系統設計和保障衛星鏈路質量至關重要。

目前國際上雨衰預測模型很多,如Crane Global模型、Assis-Einloft Improved模型、Silva Mello模型[7]、ITU-R 模型[8]、DAH模型和UK模型等。文獻[9]基于指數雨胞分布模型,提出了一種地面視距鏈路的全概率雨衰預測模型。利用ITU-R視距鏈路數據庫數據回歸得到了路徑調整因子中的參數,提高了預測精度,但對物理基礎要求較高,計算量較大;文獻[10]基于SC EXCELL模型,利用合適的測量數據作為輸入,如月降雨高度、降雨速率、降雨分布等,對預測流程作了詳盡的介紹,并將年統計數據得到的精確度與月統計數據得到的精確度進行比較,月精確度較高,但是對物理基礎要求較高,且對不同頻率雨衰變換對時間間隔要求較大,不同天氣要采取不同的采樣時間間隔,在實際應用過程中不易把握。文獻[11]基于微波雨衰的冪律關系,采用層析技術建立了2維降雨場反演模型,能夠提供高時空分辨率的2維降雨強度分布,提高了反演精度,但處理層析技術時對雨衰模型呈非線性頻段處理采用了線性化的處理方法,反演誤差較大。通過這些模型可以獲得雨衰的長期統計特性,但降雨及其產生的雨衰是不平穩的,很難精確預測短期(如下一分鐘)降雨衰減,不利于實時要求較強的業務。

在此基礎上,本文采用時間序列分析法,對降雨產生的雨衰時間序列間的相關性進行研究,提出一種基于非平穩時序ARIMA模型的W 頻段雨衰實時預測方法。本文選取羅馬地區降雨數據進行驗證,將ARIMA模型與在W頻段經過驗證的ITU-R和Silva Mello雨衰預測模型進行比較來驗證該方法的有效性。

2 時間序列分析法和ARIMA模型

上述傳統雨衰預測模型均基于對各類參數的計算基礎上,并未考慮數據之間的內部關聯性,且絕大多數模型僅驗證到55 GHz,在W頻段進行有效雨衰預測經過驗證的模型目前僅為少數,如ITU-R模型和Silva Mello模型。本文采用時間序列分析法對W頻段雨衰進行預測,并與ITU-R, Silva Mello模型預測結果進行比較。

2.1時間序列分析法

時間序列分析法[12]是利用按時間順序排列的一組數列,應用數理統計方法對序列內部的相互關系進行處理,以預測事物未來發展情況。其中,平穩時間序列的均值、方差、協方差等統計特性并不隨時間推移而變化,且這些數字特征具有遍歷性(能夠保證樣本統計量是總體參數的一致估計),在不同時點上呈現出相似的規律性,便于根據序列在過去的特征規律,建立擬合模型,對未來趨勢加以預測。

而在實際降雨衰減過程中,很多序列是不平穩的[13],其均值、方差和自協方差等數字特征會隨時間推移而變化,也就是在不同時點上的統計規律不同,且不具有遍歷性,所以不能直接利用此類序列的以往信息來預測未來的可能情景,需要采用一些方法將非平穩序列變換為平穩序列。而差分可以從緩慢的長期變化中,分離出微小的、快速的變化因素,可以將非平穩序列處理為平穩序列,特別適合處理此類問題[14]。因此本文采用ARIMA模型來研究雨衰時間序列間的相關性,以此預測未來雨衰情況。

2.2 ARIMA模型

自回歸移動平均(Auto Regression Moving Average, ARMA)模型[15]廣泛應用于物理科學、經濟學、生物學等領域隨機現象的預測,非常適合于預測平穩時間序列。假設是一組有效的平穩時間序列,則的ARMA(,)模型可表示為

ARIMA模型[15]是對ARMA模型進行優化得來的,屬于線性模型,適合于預測非平穩時間序列。在實際降雨中,絕大多數雨衰時間序列是不平穩的,在預測前需采用一些方法將非平穩序列變換為平穩序列,如差分處理。假設是一組非平穩時間序列,則的ARIMA(,,)模型可表示為

圖1為構建ARIMA模型的流程,包括以下幾步:平穩性檢驗、模型識別、參數估計、診斷驗證和預測分析。本文結合羅馬地區降雨數據,進行建模,并將該ARIMA(1,1,6)模型在不同極化方式、預測間隔和時序個數的條件下進行比較,然后分別與ITU-R, Silva Mello模型進行比較。

圖1 ARIMA模型構建流程

3 ARIMA模型構建

3.1平穩性檢驗

本文采用羅馬地區降雨數據。使用文獻[16]中的方法及SPSS(Statistical Product Service Solutions)統計軟件得到羅馬地區71.9~76.8 GHz處雨衰仿真時間序列圖,如圖2(a)所示。

采用自相關函數ACF來驗證時間序列的平穩性。統計學家給出了許多關于自相關函數的估計方法,并得出結論,延遲的自相關函數最令人滿意的估計是

采用SPSS軟件求解該序列的自相關函數ACF。若其自相關函數ACF迅速減小到零,則滿足平穩性條件。圖2(b)為使用SPSS軟件得到的初始雨衰時間序列的ACF。

由圖2(b)可看出,初始雨衰時間序列的ACF沒有迅速減小為零,說明該序列非平穩,需對其進行差分處理,使其平穩化。圖3(a)和圖3(b)分別為一階差分后的雨衰時間序列及其ACF。

由圖3(b)可知,經一階差分的雨衰時間序列的ACF迅速減小為零,說明該序列被平穩化,且差分階數=1。

圖2 初始雨衰時間序列及其ACF

圖3 差分雨衰時間序列及其ACF

3.2模型識別

檢驗平穩性和確定差分階數之后,使用自相關函數ACF和偏自相關函數PACF確定和的值,即確定ARIMA(,,)模型結構。圖4為一階差分后的平穩雨衰時間序列的抽象模型示意圖。

圖4 平穩雨衰時間序列抽象模型示意圖

圖5 差分雨衰時間序列的PACF

3.3參數估計

在3.1節中對平穩性進行了檢驗,從檢驗結果可以看出式(5)能夠使用。要想得到的估計值,就得先求出和的值。可使用最小均方估計、最小方差法、點估計或使用ACF直接估計ARIMA模型的參數。本文采用最小方差法進行計算,即選使均方偏差最小。

表1 ARIMA(1,1,6)的參數和

表1 ARIMA(1,1,6)的參數和

參數日期 2010年3月14日2010年3月29日2010年4月14日2010年5月13日 -0.8539-0.2116-0.9939 0.7422 -0.5715-0.4952 0.1541 1.0027 0.3768 0.3945-0.7961 0.0884 0.1358 0.2459 0.3441-0.0837 -0.0230 0.0059 0.1900-0.0686 -0.0613 0.0124 0.0641 0.3431 0.0191-0.0110 0.1349 0.1748

3.4診斷檢驗

可通過預測時間序列與真實序列間殘差序列的ACF和PACF,來確認該模型是否合適。

求解上述預測差分時間序列與真實序列間殘差序列的ACF和PACF,可知,ACF和PACF迅速減小為零,表明該殘差序列與白噪聲序列相似,可確認ARIMA(1,1,6)模型對羅馬地區該樣本中雨衰測量時間序列是可行的。并得出以下結論:初始雨衰時間序列不平穩;經差分的雨衰時間序列在差分階數=1時是平穩的;=1,=6可行。也就是說,ARIMA(1,1,6)模型可用于預測W頻段雨衰時間序列。得到的ARIMA(1,1,6)模型為

3.5 預測分析

誤差是不確定性最直觀的表現形式,預測結果的不確定性可由預測誤差的變化趨勢很好地展現出來。該模型通過控制不同的極化方式、預測間隔和時序個數,與ITU-R模型的預測精度進行比較,驗證本方法的有效性、參數計算簡單和高預測精度。

4 仿真分析

4.1 不同極化方式下的預測精度

從圖6(a)可以看出,分別采用水平極化、圓極化和垂直極化這幾種不同的方式進行預測,三者得到的預測精度之間的差別并不是很明顯,與此同時隨頻率的增大,三者的預測精度會逐漸朝同一程度趨近。由圖6(a)可知,在71.0~76.6 GHz范圍內,垂直極化和圓極化方式下預測精度差別較小,在76.6 ~86.0 GHz范圍內,水平極化和垂直極化方式下的誤差大致相同,但在整個W頻段范圍內,三者間預測精度差異很小,所以基本上可以忽略極化方式對該模型參數的影響,同時也可簡化該模型在不同極化方式下的預測計算量,降低了復雜度。

4.2不同預測間隔下的預測精度

選擇一階差分ARIMA(1,1,6)模型,采用垂直極化方式,時間序列個數,圖6(b)為仿真結果。由圖6(b)可以看出,不同預測間隔下,ARIMA(1,1,6)模型的預測值間的誤差隨頻率增大而升高,該結論與之前對平穩性的驗證中得出的結論是相符合的,這說明所采用的算法及進行的仿真都是可取的;頻率預測間隔不同則預測模型的預測精度就會不同,因為降低頻率預測間隔,將會增強前導雨衰數據之間的相關性,而差分是線性變換,線性變換提高了數據之間的相互聯系,所以將增強時間序列中各元素之間的相關性,從而可以提高預測精度,更精確地估計模型參數;而且預測某一頻段范圍的降雨衰減值時,頻率預測間隔大小與預測數之間的關系比較密切,預測數的多少由頻率預測間隔的大小來決定,從而導致會生成不同的累積誤差,且隨著預測間隔的增大生成的累積誤差會減小,因此當預測間隔為0.20 GHz時曲線的變化趨勢比較顯著,在71.0~75.2 GHz范圍內,預測間隔為0.10 GHz與0.05 GHz時兩者間預測精度差別較小;在75.2~86.0 GHz范圍內,預測間隔為0.20 GHz與0.05 GHz時兩者間預測精度的差別較小,同時隨頻率的增大兩者的預測精度會逐漸朝同一程度逼近。

表2模型參數

時序個數n預測間隔(GHz)偏相關系數 200.20 0.6539-0.0347-0.00780.02700.0463 0.10 0.0621 0.0717 0.08140.07260.1855 500.05-2.6132-1.5761 1.70032.54640.9676 0.10 0.8236 0.0427 0.05690.06310.1849 0.7923 0.0376 0.04850.06260.1763 0.8932 0.0016 800.10 0.0050 0.0057 0.00640.00730.0078 0.6835 0.0267 0.01570.00810.0052 0.8478 0.0593 0.06500.15790.0796 0.0083 0.0086 0.00890.07020.0815

圖6 不同極化方式和預測間隔的預測誤差比較

4.3不同時序個數下的預測精度

選擇一階差分ARIMA(1,1,6)模型,采用垂直極化方式,頻率預測間隔為0.10 GHz,圖7(a)為仿真結果。由圖7(a)可以看出,時間序列個數不同則預測模型的預測精度也會有所差異,原因是前導雨衰數據量的多少決定了時間序列個數,前導雨衰數據描述了降雨數據和雨衰數據的部分特性,它的數據越多,對數據的描述就會越充分,就越能準確地估計數據間的相關程度,從而越能精確地估計模型參數;同時注意到待計算模型參數的多少由時間序列個數的多少來確定,從而導致產生不同的累積誤差,累積誤差會隨著時間序列數目的增多而增大,所以時間序列個數為80時誤差曲線的趨勢變化較為明顯,在71~75 GHz內,時序個數為80與50時兩者預測誤差曲線基本吻合;在75~79 GHz內,時序個數為50和20時兩者間預測精度的差別較小;在79~86 GHz內,時序個數為80與20時兩者間預測精度較為接近,且誤差較小,與此同時隨著預測頻率的增大兩者的預測精度會逐漸朝同一程度逼近,且在整個W頻段范圍內,時序個數為50時預測誤差變化保持基本平穩。

4.4不同預測模型間的預測精度比較

圖7(b)為分別利用ITU-R模型、ARIMA(1,1,6)模型和Silva Mello模型對羅馬地區某一時段雨衰預測的雨衰-頻率變化曲線圖,采用垂直極化方式,其中ARIMA(1,1,6)模型采用的預測間隔為0.10 GHz,時序個數為50,從圖7(b)可看出,通過ARIMA(1,1,6)模型得到的W頻段雨衰值與ITU-R, Silva Mello模型誤差不超過1.5 dB。圖8為分別使用這3種模型得到的均方根誤差-頻率曲線圖,由圖8可看出,在W頻段,ARIMA(1,1,6)模型預測精度較高,且趨勢較吻合,具有較好的預測效果。

4.5預測序列與仿真序列吻合度分析

圖9為在76 GHz處使用ARIMA(1,1,6)模型,對2010年5月26日羅馬地區某場雨仿真時序中的前1~440個序列預測第441~590時序得到的雨衰-序列圖,其中,ARIMA(1,1,6)模型采用的預測間隔為0.10 GHz,時序個數為50。由圖9可看出,在W頻段,441~590預測序列和仿真序列間重合度較好,且變化趨勢基本吻合,說明ARIMA(1,1,6)模型預測精度較高,有較好的預測效果。

圖7 不同時間序列個數的預測誤差及不同雨衰模型預測值比較

圖8 均方根誤差與頻率的關系 圖9 2010年5月26日雨衰預測與仿真時序對照

5 結論

本文針對降雨及其產生的雨衰多為不平穩事件,提出了基于非平穩時間序列ARIMA的W頻段雨衰實時預測模型。該模型利用差分變換將非平穩時間序列轉化為平穩時間序列,將該ARIMA(1,1,6)模型在不同極化方式、預測間隔、時序個數的條件下進行比較,結果表明,極化方式對模型參數的影響基本可忽略,而在滿足預測間隔0.10 GHz,時序個數50,一階差分的條件下預測誤差不超過,然后與ITU-R, Silva Mello模型在垂直極化方式、預測間隔0.10 GHz,時序個數50的條件下進行比較,ARIMA模型屬于線性模型,適合于短期預測,而ITU-R, Silva Mello模型為長期統計模型,該3種模型雨衰隨頻率變化趨勢較吻合,而雨衰值相差不大,說明該ARIMA(1,1,6)模型與長期統計模型雨衰分布基本趨于一致。在使用該模型進行預測時,預測時序與仿真時序吻合度較高,驗證了該方法可用于W頻段雨衰預測,且預測精度較高,模型表述簡單,因此,對W頻段自適應抗衰減技術的研究和應用具有一定的參考價值。

[1] Cianca E, Rossi T, Yahalom A,. EHF for satellite communications: the new broadband frontier[J]., 2011, 99(11): 1858-1881.

[2] Nessel J A, Acosta R J, and Miranda F A. Preliminary experiments for the assessment of V/W-band links for space-earth communications[C]. IEEE Antennas and Propagation Society International Symposium(APSURSI), Orlando, FL, 2013: 1616-1617.

[3] Stallo C, Cianca E, Mukherjee S,.. UWB for multi-gigabit/s communications beyond 60 GHz[J]., 2013, 52(1): 161-181.

[4] Stallo C, Mukherjee S, Cianca E,.. System level comparison of broadband satellite communications in Ka/Q/W bands[C]. 2012 IEEE First AESS European Conference on Satellite Telecommunications (ESTEL), Rome, 2012: 1-10.

[5] Riva C, Capsoni C, Luini L,.. The challenge of using the W band in satellite communication[J]., 2014, 32(3): 187-200.

[6] Sacchi C and Rossi T. Overview of PHY-layer Design Challenges and Viable Solutions in W-band Broadband Satellite Communications[M]. Springer Berlin Heidelberg, 2010: 3-18.

[7] Pontes M S, da Silva Mello L, Willis M J,.. Eperimental data and testing procedures for modelling of propagation effects on terrestrial radio links from C to W bands[C]. 2012 6th European Conference on Antennas and Propagation (EUCAP), Prague, 2012: 86-90.

[8] Series P. Propagation data and prediction methods required for the design of Earth-space telecommunication systems[S]. Recommendation ITU-R, 2009.

[9] 趙振維, 盧昌勝, 林樂科. 基于雨胞分布的視距鏈路雨衰減預報模型[J]. 電波科學學報, 2009, 24(4): 627-631.

Zhao Zhen-wei, Lu Chang-sheng, and Lin Le-ke. Prediction model of rain attenuation based on the EXCELL rain cell model for the terrestrail line-of-sight systems[J]., 2009, 24(4): 627-631.

[10] Capsoni C and Luini L. The SC EXCELL model for the prediction of monthly rain attenuation statistics[C]. 2013 7th European Conference on Antennas and Propagation (EUCAP), Gothenburg, 2013: 1382-1385.

[11] 姜世泰, 高太長, 劉西川, 等. 基于微波鏈路的降雨場反演方法研究[J]. 物理學報, 2013, 62(15): 154303.

Jiang Shi-tai, Gao Tai-chang, Liu Xi-chuan,..Investigation of the invesion of rainfall field based on microwave links[J]., 2013, 62(15): 154303.

[12] 王輝, 魏文博, 金勝, 等. 基于同步大地電磁時間序列依賴關系的噪聲處理[J]. 地球物理學報, 2014, 57(2): 531-545.

Wang H, Wei W B, Jin S,Removal of magnetotelluric noise based on synchronous time series relationship[J]., 2014, 57(2): 531-545.

[13] Nigam R, Nigam S, and Mittal S K. Stochastic modeling of rainfall and runoff phenomenon: a time series approach review[J]., 2014, 4(2): 81-109.

[14] Das D and Maitra A. Time series prediction of rain attenuation from rain rate measurement using synthetic storm technique for a tropical location[J].-, 2014, 68(1): 33-36.

[15] 楊明, 金晨輝, 張國雙. 截斷差分概率的上界估計與應用[J]. 電子與信息學報, 2014, 36(9): 2124-2130.

Yang M, Jin C H, and Zhang G S. Evaluation and application of the upper bound probability of the truncated differential[J].&, 2014, 36(9): 2124-2130.

[16] 弓樹宏. 電磁波在對流層中傳輸與散射若干問題研究[D]. [博士論文], 西安電子科技大學, 2008: 25-61.

Gong S H. Study on some problems for radio wave propagation and scattering through troposphere[D]. [Ph.D. dissertation], Xidian University, 2008: 25-61.

[17] Ling S, Peng L, and Zhu F. Inference for a special bilinear time-series model[J]., 2015, 36(1): 61-66.

[18] 汪榮鑫. 隨機過程[M]. 西安: 西安交通大學出版社, 2006: 132-143.

Wang R X. Random Processing[M]. Xi’an: Xi’an Jiaotong University Press, 2006: 132-143.

[19] Hansen P R and Lunde A. Estimating the persistence and the autocorrelation function of a time series that is measured with error[J]., 2014, 30(1): 60-93.

Rain Attenuation Prediction at W Band Based on Non-stationary Time-series ARIMA Model

Yang Feng Xue Bin Liu Jian

(,,’710077,)

Most existing rain attenuation prediction models are only tested at 55 GHz. There are small numbers of tested W-band rain attenuation prediction models, but these tested models have issues with high complexity and large quantities of calculations. A real-time prediction method is proposed that has a simpler structure and smaller quantity of calculations. The proposed method is based on the ARIMA model, which utilizes the relationship among the time series to establish a prediction model, conducts a stationary test on the original time series, transforms the nonstationary time series into a stationary time series by using a difference transform, and estimates the parameters of the stationary time series. This sequentially transforms the traditional nonlinear prediction into a linear prediction. First, the ARIMA (1,1,6) model is compared under the conditions of different polarizations, prediction intervals, and numbers of time series. Then the proposed model is compared with the International Telecommunication Union-R (ITU-R) and the Silva Mello rain attenuation prediction models using the conditions of vertical polarization, a prediction interval of 0.10 GHz, and a number of 50 time series. Finally, the forecasting time series with the simulant series are compared. The result shows that the prediction error between the ARIMA model, the ITU-R model and the Silva Mello model does not exceed, and that the change trend of the rain attenuation is basically the same. and the goodness of fit between the forecasting and simulant time series is good, which means that the proposed model can be applied to forecasting the rain attenuation in the W band, and that it has the advantages of simpler structure and high precision in prediction.

Rain attenuation prediction; W band; Non-stationary time-series; ARIMA; Difference

TN927

A

1009-5896(2015)10-2475-08

10.11999/JEIT150472

2015-04-27;改回日期:2015-07-20;

2015-07-27

薛斌 1085028568@qq.com

航空科學基金(20120196002)和雷達型空空導彈抗拖曳式誘餌方案研究

The National Aerospace Science Foundation of China (20120196002); The Research of Anti-TRAD in Air-to-Air Missiles with Radar Seeker

楊 峰: 男,1975年生,教授,研究方向為現代無線通信、衛星與移動通信等.

薛 斌: 男,1990年生,碩士,研究方向為衛星與移動通信.

劉 劍: 男,1978年生,副教授,研究方向為現代無線通信、陣列信號處理等.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产三区二区| 极品国产一区二区三区| 在线视频97| 精品视频福利| a色毛片免费视频| a级毛片免费网站| 第一区免费在线观看| 欧美一级片在线| 欧美区一区二区三| 久久精品亚洲专区| 精品国产美女福到在线不卡f| 成人欧美日韩| 国模极品一区二区三区| 亚洲综合久久成人AV| 国产午夜人做人免费视频中文| 成年片色大黄全免费网站久久| 亚洲一区二区日韩欧美gif| 青青极品在线| 国产91视频免费| 四虎AV麻豆| 无码免费的亚洲视频| 波多野结衣的av一区二区三区| 欧亚日韩Av| 国产91全国探花系列在线播放| 国产白浆视频| 欧美不卡视频在线观看| 免费观看国产小粉嫩喷水| 国产黄在线免费观看| 黄色福利在线| 无码精品国产VA在线观看DVD| 最近最新中文字幕免费的一页| 亚洲国产精品一区二区第一页免 | 米奇精品一区二区三区| 国产凹凸一区在线观看视频| 亚洲VA中文字幕| 亚洲二区视频| 亚洲成人免费看| 伦伦影院精品一区| 亚洲天堂首页| 国产嫖妓91东北老熟女久久一| 欧美在线伊人| 欧美天堂久久| 无码日韩精品91超碰| 日本一区二区不卡视频| 无码av免费不卡在线观看| 国产三级精品三级在线观看| 国产一级毛片在线| 亚洲激情99| 99re在线视频观看| 亚洲欧洲日韩综合色天使| 国产精品爽爽va在线无码观看| 欧美一级在线播放| 国产一区二区精品福利| 麻豆国产精品| 亚洲天堂视频在线播放| 久久免费视频6| 国产高清精品在线91| 热热久久狠狠偷偷色男同| 视频二区欧美| 亚洲欧美自拍一区| 在线中文字幕网| 亚洲日韩AV无码精品| 伊人丁香五月天久久综合| 看国产毛片| 激情无码视频在线看| 国产日韩久久久久无码精品| 国产18在线播放| 亚洲欧洲日韩久久狠狠爱| 无码'专区第一页| 欧美不卡二区| 九色最新网址| 夜夜操天天摸| 米奇精品一区二区三区| 99这里只有精品6| 国产无码制服丝袜| 欧美成人精品在线| 久久国产精品77777| 国产丰满成熟女性性满足视频| 国产免费久久精品99re不卡| 青草国产在线视频| 成人午夜免费观看| 国产精品久久久久久久伊一|