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

基于貝葉斯模型組合的隨機森林預測方法

2019-04-13 13:23:22董娜常建芳吳愛國
湖南大學學報·自然科學版 2019年2期

董娜 常建芳 吳愛國

摘 ??要:為了能夠精準可靠地估計太陽能輻照度,本文提出一種基于貝葉斯模型組合的隨機森林算法用于太陽能輻照度預測. 首先,引入K-means聚類和K折交叉驗證將氣象數據訓練集生成多個訓練子集,以增加訓練子集的多樣性并保證均勻采樣. 其次,將隨機森林作為基學習器建立集成學習預測模型,導入訓練子集并訓練各個隨機森林. 之后,依據各個隨機森林在驗證集上的預測性能,采用貝葉斯模型組合算法制定組合策略. 個體隨機森林在測試集上的預測值經過模型組合策略得到最終輸出. 最后,基于氣象實測數據建立仿真實驗,并引入其他四種預測方法進行對比仿真研究,通過實驗結果驗證了文中所提出預測方法在太陽能輻照度預測問題中的準確性和可靠性.

關鍵詞:K均值聚類;交叉驗證;隨機森林;貝葉斯模型組合;太陽能輻照度

中圖分類號:TP181 ????????????????????????????????文獻標志碼:A

Random Forest Prediction Method Based on Bayesian Model Combination

DONG Na,CHANG Jianfang,WU Aiguo

(School of Electrical and Information Engineering,Tianjin University,Tianjin 300072,China)

Abstract: To accurately and reliably estimate the solar irradiance, a random forest algorithm was proposed based on the Bayesian model combination for solar irradiance prediction. Firstly, the K-means clustering and K-fold cross validation were introduced to generate multiple training subsets so as to increase the diversity of training subsets and to ensure uniform sampling. Secondly, the random forests were defined as base learners to establish an ensemble learning prediction model,with each training subset being used to train the corresponding individual random forest. Then, according to the prediction performance of each individual random forest on the verification set, the Bayesian model combination algorithm was applied to formulate the combination strategy. The prediction values of individual random forest on the test set were fused to the final output through the model combination strategy. Finally, the proposed method was applied to solve the solar irradiance prediction problem. Simulation experiments were carried out by measured meteorological data. Other four kinds of prediction methods were also introduced to establish the contrast experiments,and the accuracy and reliability of the proposed method in the solar irradiance prediction were verified by comparison results.

Key words: K-means clustering;cross validation;random forest;Bayesian model combination;solar irradiance

太陽能在光熱領域和光電領域被廣泛應用并被視為最佳代替能源. 季節、氣候、云層密度等氣候因素引起太陽能輻射量的不確定性制約了其應用領域的發展. 高精度的預測方法一直是太陽能預測研究的熱點[1-2].

當前,太陽能輻照度的預測研究主要是使用支持向量機[3-6] (SVM)和人工神經網絡[7-9](ANN)算法. 這類學習算法難于平衡訓練集的訓練誤差和測試集的泛化誤差之間的關系[10],在訓練過程中容易出現過擬合或欠擬合的現象. 然而,在太陽能供熱系統的熱水供應量估計研究中,保證預測精度的同時預測結果的可靠性顯得更為重要[11]. 集成學習為提

高預測結果的可靠性提供了思路. 集成學習(ensemble learning,EL)[12]將多個基學習器組合在一起,常可獲得比單個基學習器更顯著的泛化性能和可靠性.

本文提出一種基于貝葉斯模型組合的隨機森林預測方法(Bayesian model combination-ensemble learning,BMC-EL)用于太陽能輻照度預測,使用隨機森林作為基學習器建立集成學習模型. 首先引入K-means聚類[13]和K折交叉驗證[14]將氣象數據訓練集劃分為多組訓練子集,以增加基學習器輸入樣本的多樣性. 其次導入訓練子集并訓練各個隨機森林. 之后,依據個體隨機森林在驗證集上的預測精度,采用貝葉斯模型組合[15]算法制定個體隨機森林的組合策略. 最后將各個隨機森林在測試集上的預測輸出依據模型組合策略得到最終太陽能輻照度預

測值.

使用美國氣象協會2013-2014年太陽能預測競賽數據[16]作為數據集,使用經典ANN、SVM 、Multikernel_SVM、K-means_RBF算法建立基于氣象數據的太陽能輻照度預測對照實驗. 實驗結果驗證了提出的算法在太陽能輻照度預測研究中的準確性和可靠性.

1 ??訓練子集多樣性處理

待組合基學習器之間的差異性比較顯著時,集成學習模型會擁有更好的性能. 故增加訓練子集的多樣性以提高基學習器輸入樣本的差異性. 基于氣象數據的太陽能輻照度預測研究中,不同天氣狀況下氣象數據呈現差異性,然而傳統隨機采樣過程會導致訓練子集中不同天氣狀況樣本分布不均勻. 針對上述問題,提出K-means聚類和K折交叉驗證方法增加訓練子集的多樣性,如圖1.(為了區別K-means聚類和K交叉驗證的下標,后文中將K折交叉驗證改為M折交叉驗證)

假設需要生成M個訓練子集{D1,D2,…,DM}. 對K-means聚類生成的簇C1進行M折交叉驗證并隨機生成M個包{b1,b2,…,bM}. 將{b2,b3,…,bM}導入訓練子集D1,將{b1,b3,…,bM}導入訓練子集D2,依次將不同的M-1個包導入對應的訓練子集,直至將{b1,b2,…,bM-1}導入訓練子集DM. 類似地,對簇{C1,C2,…,Ck}都進行M折交叉驗證,并分別將其不同的M-1個包導入訓練子集{D1,D2,…,DM}.

先聚類再交叉驗證,可以使每個訓練子集中都包含不同類型天氣對應的氣象數據,這保證了均勻采樣. 交叉驗證方法劃分訓練子集增加了訓練子集的多樣性.

2 ??隨機森林基學習器

集成學習可以通過組合策略提高預測方法的可靠性. 隨機森林中回歸樹的剪枝操作可以有效降低過擬合的風險,它簡單高效,容易實現,計算開銷小,在很多分類回歸問題中展現出強大的性能. 故本文采用隨機森林算法作為基學習器.

本文采用CART回歸樹建立隨機森林的基學習器. 訓練集D={(x1,y1),(x2,y2),(x3,y3),…,(xn,yn)},輸入樣本xi = (x1i,x2i,…,xZi)包含Z個屬性變量,輸出 Y = (y1,y2,…,yn)為連續值. 回歸樹的節點對樣本xi(1 < i < n)的屬性變量j設置切分點s,該輸入變量大于s劃分為一個區域,否則劃分到另一個區域. 對劃分得到的區域使用不同的屬性變量進一步劃分,依據節點的切分點將輸入劃分為m個區域,分別記為R1,R2,…,Rm. 定義每個區域的輸出值分別為:c1,c2,…,cm. 則CART的模型為公式(1):

由上式可得,優化區域Rm的輸出值cm可以使得平方誤差最小化. 易得當cm 為屬于Rm區域的輸入樣本對應真實輸出值的均值時,平方誤差E最優,即?觬m = ave(yi|xi∈Rm).

假設選擇樣本中的變量x(j)為切分變量,節點取值s為切分點. 輸入樣本中變量j與切分點s比較就可以得到區域R1(j,s)={x|x(j)≤s}和區域R2(j,s)={x|x(j)>s}. 當j和s設為確定值時,區域R1(j,s)和R2(j,s)包含的樣本也確定. 故需要確定每個區域的輸出值c1和c2使各自區間上的平方差最小,如式(3):

則1 = ave(yi|xi∈R1(j,s)),?觬2 = ave(yi|xi∈R2(j,s)).然后遍歷樣本中所有的變量,不同切分變量的最優切分點s得到的平方誤差最小時記為最優切分變量j. 類似的,對切分好的區域進一步劃分,求取最優的切分變量和切分點,最終得到回歸樹f(x) = ■cm I(x∈Rm).

3 ??貝葉斯模型組合

3.1 ??貝葉斯模型平均

貝葉斯模型平均(Bayesian model averaging,BMA)是為解決模型的不確定性而提出的. 它是通過模型在驗證集上預測精度的后驗概率作為模型的權重,對多個隨機森林模型賦以合理的權重,解決單個模型的不確定性和單一性,將多個模型組合到一起的降低風險的方法.

給定數據集D,樣本di是由基學習器隨機森林的輸出值xi和太陽能輻照度真實值yi組成. 模型空間H是由有限個個體假設近似,h作為模型空間的個體假設. 在模型空間和數據集D條件下yi的后驗分布為:

式中:p(yi|xi,D,H)為所有個體假設估計yi的后驗分布加權平均值, 其中,p(yi|xi,h) = p(yi|θk,h,D)×

p(θk|h,D)dθk為假設空間h對yi的預測分布, θk是

個體假設h對應的參數向量.

通過BMA,數據集D下個體假設h的后驗概率(h假設作為數據生成模型的后驗概率)p(h|D)可以由式(5)計算:

p(h).p(D|h)=p(D|θk,h)p(θk|h)dθk是個體假設h的積分似然估計,p(θk|h)是h對應的向量參數θk的先驗分布,p(D|θk,h)是似然估計. p(h)是個體假設h的先驗概率. 雖然集成學習方法中引入訓練集采樣擾動和屬性擾動增加基學習器的差異性,但是為保證所有基學習器都有較高的預測性能,基學習器的初始參數設置并無差異,故先驗概率p(h)無需“偏袒”某一個個體假設,本文中p(h)=(k為假設空間中個體假設的數量).

3.2 ??貝葉斯模型組合

貝葉斯方法在理論上是最優的,并且在許多任務中具有很好的性能. 貝葉斯模型平均也被視為集成學習中結合基學習器的一種標準方法. 然而在貝葉斯模型平均中,積分似然估計的計算方式容易使輕微精度提升的假設獲得極高的權重[16],貝葉斯模型平均比stacking更容易過擬合,對模型的近似誤差非常敏感,且表現性能更差[17].

為了在太陽能輻照度預測試驗中更加高效地獲得集成學習的固有優勢,組合策略應該側重地反映各個假設空間的優勢互補,而不僅僅是通過貝葉斯模型平均找出最優的假設.

針對上述問題,為貝葉斯模型平均增加假設空間E建立貝葉斯模型組合,將公式(4)修改為式(6):

p(yi|xi,D,H,E) =p(yi|xi,H,e)p(e|D) ???(6)

式中:e是組合模型空間E中的個體假設模型. 貝葉斯模型平均和貝葉斯模型組合示意圖如圖2和圖3所示.

這樣的修正克服了貝葉斯模型平均給個體假設h所有權重的傾向.

4 ???基于貝葉斯模型組合的隨機森林預測

方法

4.1 ??預測方法流程

基于貝葉斯模型組合的隨機森林預測方法流程圖如圖4所示.

基于貝葉斯模型組合的隨機森林預測方法在太陽能輻照度預測實驗中的具體實施步驟如下.

1)首先采用公式x* = 將原始數據歸一化處理,對訓練集進行K-means聚類操作生成簇劃分{c1,c2,…,ck}. 對每個簇ck進行M折交叉驗證,并依次生成k個訓練子集{D1,D2,…,Dk}. 采用K-means聚類和M折交叉驗證劃分訓練子集,同時保證了訓練子集的多樣性和均勻采樣.

2)使用多個CART回歸樹構建隨機森林,將k個訓練子集訓練k個隨機森林算法.

3)向訓練好的k組隨機森林導入驗證集,輸出k組預測輸出值(y1,y2,…,yk). 假設驗證集的真實輸出值為y,構建矩陣(y,y1,y2,…,yk)并導入貝葉斯模型組合方法,根據k組隨機森林在驗證集的預測性能輸出模型組合策略.

4)向訓練好的k組隨機森林導入測試集,個體隨機森林輸出各自模型預測值(y1,y2,…,yk),則集成學習方法的預測輸出為p(Y|yk,D,H,E)=p(yi|xi,H,e)p(e|D).

4.2 ??預測方法的復雜度分析

CART回歸樹在尋找切分節點時需要遍歷當前特征的所有可能取值. 設數據樣本具有F個特征,每個特征有N個切分點,CART回歸樹共生成

S個內部節點,則CART回歸樹的時間復雜度為

O(F*N*S).

設每個隨機森林基學習器中包含M個CART回歸樹,集成學習預測方法中共包含了K個基學習器,則集成學習的時間復雜度為O(F*N*S*K*M).

在訓練子集采樣部分,設k個聚類中心,每個樣本包含F個特征,聚類中心的迭代次數為t,則K-means聚類的時間復雜度為O(k*F*t). 設每個簇包含m個樣本,則訓練子集采樣過程的時間復雜度為O(k*F*t+m*M).

貝葉斯模型平均包含h個體假設,驗證集包

含c個樣本,則貝葉斯模型平均的算法復雜度為

O(c*h). 貝葉斯模型組合的假設空間E包含e個體假設,則貝葉斯模型組合方法的時間復雜度為

O(c*h*e).

由于切分點個數N較大,故集成學習的時間復雜度遠大于訓練子集采樣和貝葉斯模型組合. 所以基于貝葉斯模型組合的隨機森林預測方法的時間復雜度大于集成學習方法,訓練子集采樣和貝葉斯模型組合的時間復雜度相對較小.

5 ??太陽能輻照度預測實驗

5.1 ??性能指標

均方誤差(MSE)和絕對平均誤差(MAE)作為太陽能輻照度預測的誤差評價指標,本文額外定義了平均誤差率(Average Error Rate,AER)和誤差率小于0.1的預測成功率(Rate of success,RS)兩個評價指標,如公式(7)~(9):

式中:Ypre是預測輸出;Yreal是真實值;Er是每個預測樣本的誤差率;AER為平均誤差率;Num為預測結果的樣本數;RS表示精確預測樣本的百分比,它反映了預測結果的可靠性.

5.2 ??訓練子集多樣性

將原始氣象數據歸一化處理,然后對訓練集進行K-means聚類操作. 太陽能輻照度預測實驗中將訓練集分為10個簇{C1,C2,…,C10}. 取輸入樣本的dswrf_sfc,dswrf_sfc,tmp_sfc三個屬性建立三維坐標系,簇中的樣本在坐標系中分布如圖5所示. 屬性參數范圍經過歸一化處理到[-1,1]區間.

分別對上述劃分的簇進行10折交叉驗證,將其中不同的9個包分別導入各個訓練子集. 4個訓練子集的氣象樣本分布如圖6所示. 由于聚類后又采用10折交叉驗證,訓練子集的樣本量為訓練集的90%. 由4個訓練子集的樣本分布圖可得,采樣過程并未影響不同天氣狀況的樣本分布. K-means聚類和M折交叉驗證結合的采樣過程不僅保證了均勻采樣,同時增加了訓練子集的差異性.

5.3 ??模型誤差估計及參數設置

隨機森林回歸模型是通過袋外數據(OOB)來估計模型誤差的. 隨機森林回歸模型中Bagging采樣理想狀態下會有37%的數據未被抽取,則將這些樣本進行模型的誤差估計. 由于隨機森林算法本身的屬性擾動,只有當CART回歸樹的數量達到一定量級時,隨機森林才會收斂到更低的泛化誤差. 將所有氣象預報樣本導入單個隨機森林算法,記錄OOB誤差與回歸樹數目之間的關系如圖7所示.

由圖7可知,當CART回歸樹的數量到達200時,隨機森林的OOB誤差趨于收斂. 將太陽能輻照度預測實驗中隨機森林模型CART回歸樹數目設置為200. 在太陽能輻照度仿真試驗中,集成學習模型設置10組隨機森林,對應上述10組訓練子集,對照實驗參數設置如表1所示. 對照實驗中EL預測方法訓練子集使用完整訓練集,模型組合策略采用貝葉斯模型平均,其他設置與貝葉斯模型組合的隨機森林預測方法一致.

5.4 ??實驗結果

本章將HOBA中尺度站的太陽能輻照度和周圍GEFS站點的氣象數據作為數據集(2008年之后的太陽能輻照度未公開),1994年1月1日~2004年12月31日的樣本作為訓練集. 使用隨機函數(randvector,MATLAB)為聚類和交叉驗證處理過的訓練子集重新排序,緩解相似的氣象預測樣本在訓練時接連出現. 將2005年1月1日~2006年12月31日的樣本作為驗證集. BMC-EL基于驗證集制定模型組合策略,對照實驗中利用驗證集優化預測模型的超參數. 將2007年全年的氣象預報樣本作為測試集,用于太陽能輻照度預測實驗. 對照實驗中各類預測方法的太陽能輻照度輸出值如圖8所示.

圖8中,對照實驗中各類預測方法的預測曲線和真實曲線有相近的變化趨勢. 然而各類預測方法在太陽能輻照度連續波動較大的樣本出現較大的誤差. SVM和Multikernel-SVM算法在太陽能輻照度預測中分別出現了過擬合或欠擬合. 其中,BMC-EL雖然也出現了較明顯的偏差,但由于貝葉斯模型組合策略,其預測曲線更接近真實曲線.

為了直觀地展示各類算法的預測誤差,對照實驗中太陽能輻照度預測輸出的誤差曲線如圖9所示,其中BMC-EL方法預測誤差曲線波動最小. 各類預測方法在2007年7月30日的太陽能輻照度預測時出現極大的偏差,最大的預測偏差近20 MJ·m-2,然而BMC-EL此處的偏差略小于10 MJ·m-2. 由于貝葉斯模型組合是從不同的模型空間中選擇最好的模型組合策略,故BMC-EL算法在一些較難預測的樣本上依然保證了穩定的預測精度. 貝葉斯模型組合策略極大地提高了預測方法的可靠性.

在坐標系中繪制y = x的直線表示預測輻照度和真實輻照度完全相同. 在散點圖中樣本點距離y = x直線的距離越遠則誤差越大. 將各類預測方法在測試集上的預測輸出和真實輸出繪制到坐標系中,如圖10所示.

圖10中BMC-EL預測方法對應的散點更加集中,并且更加貼近于y=x直線. 在BMC-EL的預測樣本中,在太陽能較為豐富的氣象條件下(晴,無云)散點最為集中,預測精度最高;太陽能輻照度匱乏的天氣(陰,雨,多云,儀器無讀數)包含了復雜的非線性,各類預測方法在氣象條件糟糕的樣本中都出現較大的偏差,而BMC-EL在這類樣本點的預測偏差最小,它通過貝葉斯模型組合策略有效地提高了預測方法的可靠性.

各類預測方法在測試集在不同月份的平均性能指標如表2所示.

綜述太陽能輻照度預測實驗結果,基于貝葉斯模型組合的隨機森林預測方法在太陽能輻照度預測研究中具有非常好的預測性能,可靠性強,對不同的天氣狀態下的太陽能輻照度都能實現精確可靠的預測.

6 ??總 ??結

本文提出一種基于貝葉斯模型組合的隨機森林方法用于太陽能輻照度預測,首先引入K-means聚類和M折交叉驗證將氣象數據訓練集生成多組彼此相交且不相同的訓練子集,以增加隨機森林輸入樣本的多樣性. 其次將多組訓練子集導入并訓練集成學習模型中的個體隨機森林. 之后,將多組隨機森林在驗證集上的輸出結果輸入貝葉斯模型組合算法,依據驗證集上預測性能的后驗分布制定隨機森林模型的組合策略. 最后各個隨機森林在測試集上的預測輸出經模型組合策略輸出太陽能輻照度預測值. 在太陽能仿真實驗中,BMC-EL方法通過增加貝葉斯模型組合方法顯著減少了單個隨機森林算法的不確定性,增加了太陽能預測輸出的可靠性. 多組輻照度預測實驗結果證明了所提出的預測方法預測精度高,可靠性強,可以精確可靠地預測不同氣象環境中的太陽能輻照度.

參考文獻

[1] ???田翠霞,黃敏,朱啟兵. 基于EMD-LMD-LSSVM聯合模型的逐時太陽輻照度預測[J]. 太陽能學報,2018,39(2):504—512.

TIAN C X,HUANG M,ZHU Q B. Hourly solar irradiance forecast based on EMD-LMD-LSSVM joint model[J]. Acta Energiae Solaris Sinica,2018,39(2):504—512.(In Chinese)

[2] ???路志英,任一墨,葛路琨. 基于樣條估計分位數回歸的光伏功率回歸模型[J]. 湖南大學學報(自然科學版),2017,44(10):91—98.

LU Z Y,REN Y M,GE L K. Photovoltaic power regression model based on spline estimation and quantile regression[J]. Journal of Hunan University (Natural Sciences),2017,44(10):91—98.(In Chinese)

[3] ???YANG X,JIANG F,LIU H. Short-term solar radiation prediction based on SVM with similar data[C]// Renewable Power Generation Conference. IET,2014:1.11—1.11.

[4] ???GUO W,MINGJIA L I,TAO L I,et al. Parameter identification of Hammerstein ARMAX model based on APSO-WLSSVM algorithm[J]. China Science Paper,2018,13(2):136—142.

[5] ???ALAM S,KANG M,PYUN J Y,et al. Performance of classification based on PCA,linear SVM,and multi-kernel SVM[C]//Eighth International Conference on Ubiquitous and Future Networks. IEEE,2016:987—989.

[6] ???ZHOU Y,CUI X,HU Q,et al. Improved multi-kernel SVM for multi-modal and imbalanced dialogue act classification[C]//International Joint Conference on Neural Networks. IEEE,2015:1—8.

[7] ???RABBI K M,NANDI I,SALEH A S,et al. Prediction of solar irradiation in Bangladesh using artificial neural network (ANN) and data mapping using GIS technology[C]//2016 4th International Conference on the Development in the in Renewable Energy Technology(ICDRET). IEEE,2016:1—6.

[8] ???ANAMIKA,KUMAR N,AKELLA A K. Prediction and efficiency evaluation of solar energy resources by using mixed ANN and DEA approaches[C]// Pes General Meeting | Conference & Exposition. IEEE,2014:1—5.

[9] ???YADAV A K,MALIK H,CHANDEL S S. ANN based prediction of daily global solar radiation for photovoltaics applications[C]//India Conference. IEEE,2016:1—5.

[10] ?CHEN L G,CHIANG H D,DONG N,et al. Group-based chaos genetic algorithm and non-linear ensemble of neural networks for short-term load forecasting[J]. Iet Generation Transmission & Distribution,2016,10(6):1440—1447.

[11] ?BAILI H,LI Y F. Online reliability prediction of energy systems with wind generation[C]//International Midwest Symposium on Circuits and Systems. IEEE,2016:1—4.

[12] KROGH A,VEDELSBY J. Neural network ensembles,cross validation and active learning[C]//International Conference on Neural Information Processing Systems. MIT Press,1994:231—238.

[13] ?ALOISE D,DESHPANDE A,HANSEN P,et al. NP-hardness of Euclidean sum-of-squares clustering[J]. Machine Learning,2009,75(2):245—248.

[14] ?劉芳,夏洪山,艾軍,等. 基于氧化動態模型的瀝青熱氧老化性能預測[J]. 湖南大學學報(自然科學版),2018,45(1):136—141.

LIU F,XIA H S,AI J,et al. Prediction of asphalt thermal oxidative aging performance based on oxidation dynamic model[J]. Journal of Hunan University (Natural Sciences),2018,45(1): 136—141. (In Chinese)

[15] ?MONTEIH K,CARROLL J L,SEPPI K,et al. Turning Bayesian model averaging into Bayesian model combination[C]//International Joint Conference on Neural Networks. IEEE,2011:2657—2663.

[16] AMS 2013-2014 Solar energy prediction contest,forecast daily solar energy with an ensemble of weather models[EB/OL]. https://www.kaggle.com/c/ams-2014-solar-energy-prediction-contest.

[17] ?CLARKE B. Comparing Bayes model averaging and stacking when model approximation error cannot be ignored[J]. Journal of Machine Learning Research,2003,4(4):683—712.

主站蜘蛛池模板: 成人无码一区二区三区视频在线观看| 国产精品人人做人人爽人人添| 国产麻豆aⅴ精品无码| 亚洲国产精品日韩av专区| 日本欧美一二三区色视频| 五月天久久综合国产一区二区| 狠狠做深爱婷婷久久一区| 亚洲午夜福利精品无码不卡| 国产激情无码一区二区APP| 国产精品成人免费视频99| 亚洲综合色婷婷| 无码日韩人妻精品久久蜜桃| 亚洲天堂网在线播放| 日本成人精品视频| 日本黄色不卡视频| 精品撒尿视频一区二区三区| 伊人大杳蕉中文无码| 日韩欧美中文亚洲高清在线| 丰满人妻被猛烈进入无码| 91在线播放免费不卡无毒| 国产精品综合色区在线观看| 久久国产精品无码hdav| 亚洲婷婷丁香| 高h视频在线| 亚洲丝袜中文字幕| 日韩色图在线观看| 国产成人永久免费视频| 欧美人人干| 国产精品成人免费视频99| 97视频在线观看免费视频| 国产精品视频免费网站| 女人一级毛片| 日日拍夜夜嗷嗷叫国产| 免费一级无码在线网站| 真人高潮娇喘嗯啊在线观看| h网址在线观看| 成人午夜福利视频| 秘书高跟黑色丝袜国产91在线| 四虎成人在线视频| 国产成人一区| 97久久免费视频| 91九色视频网| 免费观看男人免费桶女人视频| 激情乱人伦| 免费va国产在线观看| 老熟妇喷水一区二区三区| 一级一毛片a级毛片| 亚洲国产成人精品无码区性色| 久久国产热| 精品久久久久久中文字幕女| 一级爱做片免费观看久久| 正在播放久久| 国产激情第一页| 久久久久久久久亚洲精品| 日本不卡在线视频| 亚洲国产精品美女| 国产高清无码第一十页在线观看| 欧美特黄一级大黄录像| 美女内射视频WWW网站午夜| 大陆国产精品视频| 男女男免费视频网站国产| 特级毛片免费视频| 亚洲欧美成人影院| 精品久久久久无码| 国产乱子伦无码精品小说| 欧美福利在线| 欧美日韩一区二区在线免费观看| 国产尤物视频在线| 亚洲成人一区二区| 国产在线拍偷自揄观看视频网站| 中文字幕欧美日韩| 久久国产av麻豆| 日韩欧美中文在线| 免费观看男人免费桶女人视频| 麻豆AV网站免费进入| 一区二区午夜| 国产人免费人成免费视频| 在线观看免费AV网| 成人福利在线视频免费观看| 亚洲男人天堂久久| 国产极品嫩模在线观看91| 欧美亚洲一区二区三区在线|