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

對未來兩個太陽周太陽活動參數的統計預測

2016-10-10 01:49:26丁煌廖云琛肖子牛

丁煌廖云琛肖子牛,

(1 中國電力科學研究院,南京 210003; 2 南京信息工程大學,南京 210044;3 中國科學院大氣物理研究所,北京 100029)

?

對未來兩個太陽周太陽活動參數的統計預測

丁煌1廖云琛2肖子牛2,3

(1 中國電力科學研究院,南京 210003; 2 南京信息工程大學,南京 210044;3 中國科學院大氣物理研究所,北京 100029)

利用支持向量機(Support Vector Machine,SVM)和后向傳播(Back Propagation,BP)神經網絡的方法,結合前23個太陽周期(1700—2008年)的周期特征數據,對第24和第25個太陽周的各個周期特征進行了預測,并且通過交叉驗證算法得出兩種方法都可達到最優。通過分析SVM與BP神經網絡方法的預測結果,均表明第25個太陽周將會達到較強的強度,且大于第24個太陽周。此外,兩種方法都預測出第25個太陽周的太陽黑子數在谷值年維持異常偏低,周期長度都會維持在10年左右。根據第24個太陽周已經過去的特征驗證,BP神經網絡的結果與實際情況更為接近,預測太陽活動在2020開始進入第25個太陽周,在2025年達到峰值,峰值年強度比第24個太陽周偏強。

太陽黑子數,太陽周期,支持向量機,后向傳播神經網絡

0 引言

眾所周知,太陽輻射變化對地球氣候的形成和長期演變有著重要的影響[1]。IPCC報告[2]指出,太陽直接輻射對氣候變化的影響雖然很小,但太陽活動的影響有可能在氣候系統相互作用的某些環節存在調制作用,使得太陽對氣候的影響超過人們的想象。Stauning[3]利用太陽周平均的黑子數代表太陽活動的強弱,發現全球溫度的變化與太陽活動的水平有聯系。事實上,許多研究也指出太陽活動對地球氣候系統有一定的影響[4-8]。基于觀測和數值模式試驗的眾多研究分析結果也揭示出地球上有許多區域存在對太陽變化和周期信號的顯著響應[9]。Misios等[10]指出,太陽活動高值年,赤道西太平洋的深對流區向東移動;Meehl等[5]指出在太陽峰值年,整個熱帶太平洋地區會出現類冷事件,而峰值年之后一至兩年,會出現類暖事件。屠其璞等[11]在分析了太陽黑子活動與中國極端事件的關系后發現,太陽黑子活動極大值年附近我國旱、澇事件出現頻率呈現顯著增加趨勢;Zhao等[12-13]確認了東亞夏季風爆發期季風區雨帶緯度位置年代際變化在一定程度上與太陽黑子周期位相有關。1991年,Friis-Christensen等[14]的研究發現,每個太陽周的時間長度和全球地表溫度有密切的聯系。最近,Maliniemi等[15]指出,北半球的冬季氣溫異常型在太陽周的不同相位(谷值期、上升期、峰值期和下降期)有顯著的差異,尤其在下降時期有明顯的一致性特征。這些結果說明,太陽活動在一個太陽周內可以表現為周期強度、周期長度、上升年份、下降年份、峰值和谷值等不同的特征,而這些特征的變化都會對地球的氣候產生一定的影響。因此,對這些特征的預測對于預估未來地球氣候的變化具有重要的意義。對未來太陽周的預測,一直是人們關注的問題,迄今為止還沒有物理上明確的預測方法,主要以統計方法對周期的強度和峰值的出現時間進行預測[16-21]。本文利用支持向量機(Support Vector Machine,SVM)和后向傳播(Back Propagation,BP)神經網絡的方法,結合前23個太陽周期(1700—2008年)的周期特征(周期長度、周期平均黑子數、周期峰值黑子數、周期谷值黑子數、上升年數(不包含峰值年)、下降年數(不包含谷值年))數據,對第24和第25個太陽周的各個周期特征進行了預測。本文在后文中介紹了使用的方法以及資料,并且給出了預測結果。

1 資料與方法

1.1太陽黑子資料

使用太陽影響數據分析中心(SIDC)所提供的逐年平均太陽黑子數資料計算了1700—2008年共23個太陽周期(一個太陽周期定義為谷值年到下一個谷值年的時間長度(不包含后一個谷值年))的周期特性,包括:周期長度、周期平均黑子數、周期峰值黑子數、周期谷值黑子數、上升年數(不包含峰值年)、下降年數(不包含谷值年)。接下來,我們將使用SVM以及BP神經網絡來預測第24和第25個太陽周的以上周期特性。并按照周期特性來還原出第24和第25個太陽周的太陽黑子數的時間序列。

1.2支持向量機

SVM是一種建立在統計學習理論的VC維理論和結構風險最小原理基礎上的,在解決小樣本、非線性及高維模式識別中能夠表現出許多特有的優勢,并能夠推廣應用到函數擬合等其他機器學習問題中的方法[22]。其基本思想是通過用內積函數定義的非線性變換將輸入空間變換到一個高維空間,在這個高維空間中尋找輸入變量和輸出變量之間的一種非線性關系[23-26]。這樣在高維特征空間的線性回歸就對應于低維輸入空間的非線性回歸,也使得計算的復雜度不再取決于空間維數,而是取決于支持向量的個數。

根據以上理論可知支持向量機的回歸就是為了尋找一個f(x)=ω·x+b,使其結構風險最小化。有樣本數據((xi, yi),i=1, 2, 3…, n),xi∈Rd為輸入,yi∈R為對應輸出,d為輸入維數,n為樣本總數,假定存在函數f在ε精度能夠估計所有的(xi, yi)數據,尋找回歸函數可由以下最優化問題轉化:

上式中w為回歸函數法向量;ε為損失函數,表示函數擬合精度;為松弛變量,可以放寬函數優化條件;常數C>0控制對訓練誤差超出誤差帶ε樣本的懲罰程度,即為懲罰系數。通過引入拉格朗日乘子αi和αi*(0<αi,αi*

滿足Mercer條件[26]的任何對稱的核函數都對應于特征空間的點積。核函數的引入,避免了對顯式非線性映射Φ的尋找。有許多種核函數都滿足此條件,如多項式核函數、徑向基核函數和柯西核函數等。其中徑向基核函數(式3)因其優秀的局部逼近特性在SVM中應用最為廣泛,它利用局部接受域完成函數映射:只有當輸入落入輸入空間的一個局部區域時,基函數才產生一個重要的非零響應,而其他情況下的響應近似為零。文中采用核函數即為徑向基核函數。SVM回歸結構與一個三層前饋神經網絡相似:其中輸入層基本相同,隱含層節點對應于一個輸入樣本與一個支持向量的內積核函數,輸出層節點對應于隱含層輸出的線性組合。圖1為支持向量機回歸原理示意圖。

圖1 支持向量機示意圖Fig. 1 Schematic diagram of SVM

1.3后向傳播神經網絡

后向傳播算法[28]的特點是信號的前向計算與誤差的反向傳播。信號前向計算時,輸入樣本首先由輸入層傳入,經過各隱含層(一層或多層)逐層傳輸并處理后,最后傳向輸出層,由該層輸出。在此過程中若是最后由輸出層輸出的實際值與期望值不符,則轉入誤差的方向傳播階段。誤差反向傳播可以簡單理解為將誤差作為輸入傳入各層不斷調整權值以達到修正的目的。具體做法是將誤差通過輸出層,按誤差梯度下降的方式通過隱含層向輸入層逐層反傳,并調整各層權值。信號前向計算與誤差反向傳播不斷進行,既是權值不斷調整的過程亦是網絡訓練過程,調整權值的原則是使誤差不斷地減小。此過程一直持續到網絡輸出的誤差減少到初始設定參數的程度抑或到達最大訓練次數。

圖2[29]是一個典型的三層BP神經網絡,只有一個隱含層,也可以根據實際情況適當增加隱含層層數。一般使用雙極性Sigmoid函數:

作為激發函數。由圖可見BP神經網絡包括三層,從左到右分別是輸入層、隱含層和輸出層,若有樣本數據((Pi, Ti), i=1, 2, 3,…, n),則其中的Pi作為輸入向量。首先賦予各權值的初始值,即輸入層至隱含層的連接權ωij(i=1, 2,…, n;j=1, 2,…, p)、隱含層至輸出層的連接權vjt(j=1, 2 ,…, p;t=1, 2,…, q)、隱含層各單元的輸出閾值θj(j=1, 2,…, p)和輸出層各單元的輸出閾值γt(t=1, 2,…, q)賦予區間(-1, 1)內的隨機值。之后網絡隨機選取一組輸入樣本和目標樣本提供計算。有了上述樣本等值就可以計算中間隔層單元的輸入sj,它需要有輸入樣本Pk、連接權ωi和閾值θj,計算好以后將sj通過傳遞函數來算出隱含層各單元的輸出bj。其中,

接著計算輸出層各單元的輸出Lt,過程跟sj類似,需要用到隱含層的輸出bj、連接權vjt和閾值γt,然后通過傳遞函數計算輸出層各單元的實際輸出Ct。其中,

然后可以通過以上條件結合下式獲得輸出層的各單元一般化誤差:

由式(9)可知,需要通過網絡目標向量和網絡的實際輸出計算得出。此時進入了誤差反傳的階段,有了輸出層的一般化誤差,再結合連接權vjt和隱含層的輸出bj共同得到隱含層各單元的一般化誤差:

現在有了各層誤差,即可以通過誤差來修正各連接權和閾值。先根據(11)、(12)兩式,將輸出層各單元的一般化誤差與隱含層各單元的輸出bj來修正連接權vjt和閾值γt:

其中,。再利用(13)和(14)兩式來修正連接權ωij和閾值θj,此過程需要用到隱含層各單元的一般化誤差和輸入層個單元的輸入。

式中,。此時一次訓練已經完畢,接下來隨機選取下一個訓練樣本向量提供給網絡,返回向量輸入開始初始運算,反復訓練直到滿足規定誤差范圍內,即網絡收斂。若是學習次數大于最大訓練次數,網絡則無法收斂。無論網絡收斂與否,運行到此時整個網絡學習都將結束。

圖2 BP神經網絡示意圖Fig. 2 Schematic diagram of BP neural network method

1.4交叉驗證算法及過程

采用交叉驗證的思想可以在某種意義下得到最優的建模參數,可以有效地避免過學習和前學習狀態的發生,最終能夠在有效樣本的情況下得到較理想的預報值。

所謂交叉驗證的基本思想是在一定意義下,將原始樣本進行分組,一部分作為訓練樣本,另一部分作為檢驗樣本。大致方法是,首先利用訓練樣本對模型進行訓練,然后利用檢驗樣本來檢驗由訓練樣本訓練好的模型,以預報準確率作為檢驗該模型效果的指標。本文利用K次交叉驗證(K-Fold Cross Validation),將初始樣本分割成K個子樣本,選取K個子樣本中的一個單獨的子樣本用作檢驗模型,而其他(K-1)個子樣本用作訓練模型。該交叉驗證需要重復K次,被分的每個子樣本都需要檢驗一次,最后平均K次的結果得到最后結果。隨機產生的子樣本能夠被同時重復利用于訓練和檢驗是該方法的最大優勢。將每次的結果檢驗一次,也不是次數越多越好,要根據具體情況具體采用。試驗結果表明使用K次交叉驗證能夠更好地尋找到建模最優參數。

在建模選擇參數時使用到了10次交叉驗證,具體做法是在訓練預報模型前,隨機將樣本數據集切割成10個子樣本,順序抽取一個子樣本為檢驗樣本,其他9個則作為訓練樣本來訓練模型,訓練完畢后使用選取出的檢驗樣本驗證該模型,如此循環。最后綜合對10個子樣本驗證的結果,選取最優參數和最優模型,使結果更好。

1.5建模

在建模中,為了消除數據對模型的影響,所有帶入模型中的數據都需要進行歸一化處理,處理方法為x=(xi-m)/(M-m),其中xi為因子的實際值,x為歸一化后的值,M和m分別為該因子所有取值中的最大值和最小值。本文利用SVM對逐周期太陽周期特征(周期長度、周期平均黑子數、周期峰值黑子數、周期谷值黑子數、上升年數(不包含峰值年)和下降年數(不包含谷值年))數據進行預測,最優訓練數據個數為5個。

建模之前,進行了部分敏感性試驗。由于數據是一組時間序列,考慮到訓練數據個數和訓練參數的選取對模型建立影響較大,故將數據分成前3個數據為訓練、1個數據為檢驗;4個數據為訓練、1個數據為檢驗;5個數據為訓練、1個數據為檢驗,…,14個數據為訓練、1個數據為檢驗,利用交叉驗證算法得到逐周期太陽周期特征的幾個數據集,利用5個數據作為訓練都能得到最優。

2 結果

通過使用SVM與BP神經網絡方法,得到了第24和第25個太陽周的各個周期特征(周期長度、周期平均黑子數、周期峰值黑子數、周期谷值黑子數、上升年數(不包含峰值年)和下降年數(不包含谷值年))。通過交叉驗證算法發現結果都能得到最優。可見,兩種預測模型都比較有效,其預測結果分別如表1和表2所示。BP神經網絡得到的第25個太陽周的強度弱于SVM的結果。此外,兩種方法都預測出第25個周期的太陽黑子數谷值較往年低。兩種方法的周期長度預測結果表明,第24和第25個太陽周的周期長度較短,都少于11年。同時,SVM預測的第24和第25個太陽周的上升年數較短,僅僅維持3~4年。第25個太陽周雖然在峰值年太陽黑子數較多,可以達到187個,但其低谷年太陽活動異常平靜,在整個第25個太陽周中,年平均太陽黑子數僅為61.5個,太陽活動整體看仍然維持較低的水平。BP神經網絡預測結果整體類似,第24和第25個太陽周上升年數分別為5.9和5.0a,下降年數分別為6.4和5.4a。第24和第25個太陽周在峰值年太陽黑子數均較多,其中在第25個太陽周達到247個,谷年太陽活動均比較弱,太陽黑子數不到10個,在第24和第25個太陽周,年平均太陽黑子數分別為73.4和75.6個。

表1 使用SVM方法預測的第24和第25個太陽周的各個周期特征Table 1 The 24thand 25thsolar cycles' periodic characteristics predicted by SVM method

表2 使用BP神經網絡方法預測的第24和第25個太陽周的各個周期特征Table 2 The 24thand 25thsolar cycles' periodic characteristics predicted by BP neural network method

圖3 使用SVM方法預測的逐年太陽黑子數序列(虛線后為預測的第24—25周期,標有實心圓的實線為實際太陽黑子數序列)Fig. 3 The time series of sunspot number predicted by SVM method (The predicted 24thand 25thsolar cycles are after the dotted line, and the solid line with filled circle is the actual time series of sunspot number)

通過對第24和第25個太陽周預測得到的周期特征的條件限制,再通過插值可以方便地重構第24—25周的太陽黑子數的時間序列。圖3和圖4給出了1990—2015年實際太陽黑子數的變化曲線(標有實心圓的實線)以及2008年后使用兩種方法預測的第24—25周的太陽黑子數的時間序列(虛線之后)。從圖3可以看到,SVM方法預測的第24個太陽周在2012年達到峰值,并在2019年觸及太陽活動谷底,開始進入第25個太陽周。第25個太陽周預計在2022年達到峰值,其強度會超過第24個太陽周的峰值強度。而從實際的2008—2015年的太陽黑子數序列上可以看出,第24個太陽周是在2014年達到峰值,其上升年數較SVM方法預測要多,整個周期較為推后。基于BP神經網絡方法的預測中,第24個太陽周在2014年達到峰值,在2020年進入谷底且開始第25個太陽周,并預計在2025年達到峰值。BP神經網絡預測的峰值年與實際情況一致,都位于2014年。此外,SVM與BP神經網絡方法的預測結果都表明第25個太陽周未來達到峰值的強度將大于第24個太陽周,而第24個太陽周的峰值較第23個太陽周弱,這與實際太陽黑子數的序列比較一致。通過BP神經網絡得到的預測結果與Schatten[30]的結果較為一致,都表明在2020年將會開始第25個太陽周。苗娟等[21]也指出,預計第25個太陽周從2020年開始,其強度與第23個太陽周類似,比第24個太陽周強。通過計算得出2009—2015年實際太陽黑子數與2009—2015年預測(SVM、BP神經網絡)太陽黑子數的相關系數,結果如表3所示。BP神經網絡得到的預測結果與實際太陽黑子序列的相關系數更高,通過了0.01的信度檢驗。因此,從已經過去的第24個太陽周的情況來看,BP神經網絡方法的預測結果更接近實際情況。

表3 2009—2015年實際太陽黑子數(SSN)與2009—2015年預測(SVM,BP神經網絡)太陽黑子數(SSN)的相關系數Table 3 The correlation coefficients between the 2009-2015 actual sunspot number (SSN) and the 2009-2015 predicted (SVM, BP neural network) SSN

圖4 使用BP神經網絡預測的逐年太陽黑子數序列(虛線后為預測的第24—25周期,標有實心圓的實線為實際太陽黑子數序列)Fig. 4 The time series of sunspot number predicted by BP neural network method (The predicted 24thand 25thsolar cycles are after the dotted line, and the solid line with filled circle is the actual time series of sunspot number)

3 總結

本文基于SVM和BP神經網絡方法對未來第24和第25個太陽周的主要周期特征進行了預測,結果表明兩種方法都對太陽活動的演變趨勢具有一定的預報能力,其中BP神經網絡方法更為接近實際情況。未來預計在2020開始進入第25個太陽周,其上升時間較短,2025年達到峰值,且其峰值年強度比第24個太陽周要強。

[1]趙亮, 徐影, 王勁松, 等. 太陽活動對近百年氣候變化的影響研究進展. 氣象科技進展, 2011, 1(4): 37-48.

[2]IPCC. IPCC Fourth Assessment Report: Climate Change 2007(AR4). IPCC, 2007.

[3]Stauning P. Solar activity-climate relations: a different approach. J Atmos Solar-Terrestrial Phys, 2011, 73(13): 1999-2012.

[4]Haigh J D. Te impact of solar variability on climate. Science, 1996,272(5264): 981-984.

[5]Meehl G A, Arblaster J M. A lagged warm event-like response to peaks in solar forcing in the Pacific region. J Climate, 2009, 22(13): 3647-3660.

[6]Gray L J, Beer J, Geller M, et al. Solar influences on climate. Rev Geophys, 2010, 48(4): 1032-1047.

[7]Roy I, Haigh J D. Solar cycle signals in sea level pressure and sea surface temperature. Atmos Chemi and Phys, 2010, 10(6):3147-3153.

[8]Bal S, Schimanke S, Spangehl T, et al. On the robustness of the solar cycle signal in the Pacific region. Geophysical Research Letters,2011, 38(14): 130-137.

[9]肖子牛, 鐘琦, 尹志強, 等.太陽活動年代際變化對現代氣候影響的研究進展. 地球科學進展, 2013, 28(12): 1335-1348.

[10]Misios S, Schmidt H. Mechanisms involved in the amplification of the 11-yr solar cycle signal in the tropical Pacific Ocean. J Climate,2012, 25(14): 5102-5118.

[11]屠其璞. 近百年來我國降水量的變化. 南京氣象學院學報, 1987,10(2): 177-187.

[12]Zhao L, Wang J S, Zhao H J. Signature of the solar cycle on decadal variability in monsoon precipitation over China. J Meteor Soc Japan, 2012, 90: 1-9.

[13]Zhao L, Wang J S. Robust response of the East Asian monsoon rainband to solar variability. J Climate, 2014, 27(8): 3043-3051.

[14]Friis-Christensen E, Lassen K. Length of the solar cycle: an indicator of solar activity closely associated with climate. Science,1991, 254(5032): 698-700.

[15]Maliniemi V, Asikainen T, Mursula K. Spatial distribution of Northern Hemisphere winter temperatures during different phases of the solar cycle. J Geophys Res (Atmospheres), 2014, 119(16):9752-9764.

[16]王家龍. 第 22 和 23 太陽周太陽活動預測. 天體物理學報, 1992,12(4): 369-373.

[17]王家龍. 相似周方法對第23周黑子數預測的檢驗. 第十屆全國日地空間物理學術討論會, 2003.

[18]Wang J L, Gong J C, Liu S Q, et al. Te prediction of maximum amplitudes of solar cycles and the maximum amplitude of solar cycle 24. Chinese J Astro Astrophys, 2002, 2(6): 557.

[19]柳士俊, 俞小鼎, 陳永義. 太陽黑子活動的守恒量預報. 科學通報, 2003, 48(17): 1832-1835.

[20]Du Z L, Wang H N, Zhang L Y. A running average method for predicting the size and length of a solar cycle. Chinese J Astro Astrophys, 2008, 8(4): 477.

[21]苗娟, 龔建村, 李志濤, 等. 第 25 太陽周太陽黑子數峰值預測. 中國科學(物理學, 力學, 天文學), 2015 (9): 95-101.

[22]丁煌, 陶樹旺, 肖子牛, 等. 基于WRF和SVM方法的風電場功率預報技術研究. 高原氣象, 2013, 32(2): 581-587.

[23]王定成, 方廷健, 唐毅, 等. 支持向量機回歸理論與控制的綜述.模式識別與人工智能, 2003, 16(2): 192-197.

[24]馮漢中, 徐會明, 徐琳娜. SVM方法與長江上游降水落區預報.高原氣象, 2004, 23(1): 63-68.

[25]張學工. 關于統計學習理論與支持向量機. 自動化學報, 2000,26(1): 32-42.

[26]Trevor H, Robert T, Jerome F. 統計學習基. 范明,柴玉梅,昝紅英,譯. 北京: 電子工業出版社, 2003.

[27]Janjic Z, Black T, Pyle M, et al. High resolution applications of the WRF NMM. 21st Conf on Weather Analysis and Forecasting/17th Conf on Numerical Weather Prediction, 1-5 August 2005, Washington DC, Amer Meteor Soc.

[28]韓爽. 風電場功率短期預測方法研究. 北京: 華北電力大學,2008.

[29]傅薈璇, 趙紅. MATLAB神經網絡應用設計. 機械工業出版社,2010.

[30]Schatten K. Fair space weather for solar cycle 24. Geophysical Research Letters, 2005, 32(21): 365-370.

The Statistical Prediction of Sunspot Number for the Next Two Solar Cycles

Ding Huang1, Liao Yunchen2, Xiao Ziniu2,3
(1 China Electric Power Research Institute, Nanjing 210003 2 Nanjing University of Information Science and Technology, Nanjing 210044 3 Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029)

In this work, the 24thand 25thsolar cycles' periodic characteristics were predicted by using support vector machine(SVM) method and back propagation (BP) neural network method respectively, based upon the data of previous 23 solar cycles. The authors found that both of the methods reach the optimal value after applying the cross validation algorithm, it reveals that the intensity of solar activity in the 25thsolar cycle will be enhanced comparing with the 24thsolar cycle in both the SVM and the BP Neural Network prediction; that the sun spot number (SSN) will maintain low value in the valley phase; and the cycle length will be around 10 years in the 25thsolar cycle. According to the periodic characteristics of the 24thsolar cycle that was past, there sults from BP Neural Network method prediction are more close to actual situation than that from the SVM. The results indicate that solar activity will enter 25thsolar cycle in 2020, and will reach the peak in 2025; the amplitude of peak year in the 25thsolar cycle will be stronger than that in the 24thsolar cycle.

sun spot number, solar cycle, support vector machine, back propagation neural network

10.3969/j.issn.2095-1973.2016.04.003

2016年2月18日;

2016年3月23日

丁煌(1988—),Email: dinghuang@epri.sgcc.com.cn

肖子牛(1965—),Email: xiaozn@lasg.iap.ac.cn

資助信息:國家重大科學研究計劃(2012CB957804)

主站蜘蛛池模板: 91九色国产porny| 精品视频91| 国产又色又刺激高潮免费看| 亚洲无码精彩视频在线观看| 欧美一级在线播放| 久久中文字幕不卡一二区| 666精品国产精品亚洲| 亚洲欧洲日产国码无码av喷潮| 欧美激情视频一区二区三区免费| 亚洲成在线观看 | 国产丝袜丝视频在线观看| 国产精品成人不卡在线观看| 国产一级毛片yw| 亚洲永久色| 国产毛片网站| 青草免费在线观看| 亚洲AV无码久久精品色欲| 欧美性色综合网| 亚洲欧美另类日本| 国产拍在线| 呦系列视频一区二区三区| 毛片基地美国正在播放亚洲 | 精品99在线观看| 99成人在线观看| 天堂亚洲网| 久久五月视频| 91 九色视频丝袜| 国产区免费| 国产一区二区三区夜色| 九九视频免费看| 色妞永久免费视频| 国产亚洲视频中文字幕视频 | 欧美激情第一欧美在线| 国产极品美女在线| 无码免费视频| 免费看的一级毛片| 高清欧美性猛交XXXX黑人猛交| 国产高清毛片| 免费毛片视频| 国产最新无码专区在线| 一区二区三区国产精品视频| 久久99国产综合精品女同| а∨天堂一区中文字幕| 国产精品区视频中文字幕| 成人午夜亚洲影视在线观看| 伊人丁香五月天久久综合| 国产小视频a在线观看| 日本国产精品一区久久久| 国产毛片网站| 久热中文字幕在线| 中文国产成人久久精品小说| 在线高清亚洲精品二区| 亚洲aⅴ天堂| lhav亚洲精品| 热99re99首页精品亚洲五月天| 色综合久久无码网| 亚洲欧美综合另类图片小说区| 欧美在线网| 在线毛片网站| 亚洲一级毛片免费看| 99久久精品国产综合婷婷| 国内精品视频| 亚洲欧美日韩精品专区| 国产网友愉拍精品| 免费又黄又爽又猛大片午夜| 又爽又黄又无遮挡网站| 黄色网址手机国内免费在线观看| 99热这里只有免费国产精品| 国产精品所毛片视频| 国产精品密蕾丝视频| 亚洲欧美精品在线| 天堂成人在线| 婷婷色在线视频| 1024你懂的国产精品| 91久久夜色精品| 中文字幕一区二区视频| 人妻精品久久无码区| 国产成人在线无码免费视频| a毛片免费在线观看| 激情网址在线观看| 国产精品无码翘臀在线看纯欲| 成人亚洲国产|