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

BP與LSTM神經網絡在福建小流域水文預報中的應用對比

2020-03-11 05:06:30崔巍顧冉浩陳奔月王文
人民珠江 2020年2期
關鍵詞:模型

崔巍,顧冉浩,陳奔月,王文

(1.河海大學水文水資源與水利工程國家重點實驗室,江蘇南京210098; 2.福建省水文水資源勘測局閩江河口水文實驗站,福建福州350001)

水文預報模型可以粗略地分成過程驅動模型和數據驅動模型兩大類[1],過程驅動模型以水文學概念為基礎,模擬產流和河道演進過程,從而進行流量過程預報;數據驅動模型則基本不考慮水文過程的物理機制,通過建立輸入、輸出數據之間的最優數學關系,實現流量過程的預報。隨著下墊面觀測數據、水文氣象觀測數據的不斷豐富,以及計算機計算能力的提升,數據驅動模型在水文預報中得到越來越廣泛的關注和應用。

數據驅動模型中,人工神經網絡因為具有良好的非線性映射能力,在水文領域得到了廣泛應用。以Romelhart和Mcclelland[2]在1986年提出的誤差反向傳播算法(BP算法)為基礎構建的BP神經網絡是目前最常用的一種神經網絡模型。由于其結構簡單,并且具有良好的非線性和泛化逼近能力,以及其自組織、自適應和容錯性,20世紀90年代以來就被很多研究者應用于洪水預報[3-4]及中長期徑流預報[5-6]。但BP神經網絡存在著網絡結構難以確定、訓練速度慢和預測精度低等方面的問題[7]。Elman等[8]提出的循環神經網絡(RNN)可以被看作在隱含層進行了時間上的折疊,從而可以擁有無限多個隱含層,通過在隱含層的迭代計算,使新的信息被添加到每一層中,通過無限制次數的更新,循環神經網絡可以將信息傳遞下去,使得神經網絡獲得對時間序列數據狀態特征的記憶能力[9]。RNN在流域降雨徑流預報中得到一定程度的應用,包括徑流漲退過程預測[10]、降雨徑流日尺度模擬[11]以及洪水預報[12]等。但Bengio等[13]發現RNN存在梯度消失問題(即梯度下降過程中,梯度隨著神經網絡層數的加深快速衰減,迅速接近于0,導致權重無法更新)與梯度爆炸問題(即由于初始化權重過大,誤差梯度在權重更新中累積,導致網絡權重大幅更新,權重出現快速增長,使得網絡變得不穩定),使得循環神經網絡很難學習到長周期序列數據的狀態特征。為了解決這一問題,Hochreiter和Schmidhuber[14]在一般RNN神經網絡的基礎上,提出了LSTM(長短期記憶)神經網絡,將記憶單元添加到RNN隱藏層的神經單元中,從而使得時間序列上的記憶信息變得可以控制,能夠有效解決信息的長期依賴,避免梯度消失或爆炸。近幾年來LSTM神經網絡成為一個非常受重視的流域降雨徑流預報方法。Shi等[15]將LSTM模型運用于降雨臨近預報,表明LSTM模型能夠很好地捕捉時空相關性;馮鈞等[16]則將LSTM模型和BP神經網絡模型進行出來對比應用,并將二者結合起來形成組合預報模型,進行12 h預見期的洪水預報;Tian等[17]比較了LSTM、ERNN(Elman循環神經網絡)、ESN(Echo State Network回聲狀態網絡)、NARX(動態時間序列神經網絡)和GR4J模型的日尺度徑流預報結果,發現LSTM神經網絡在面積較小的流域中有更好的預報精度;Kratzert等[18]將LSTM模型應用于受融雪影響的流域,說明LSTM相比于SAC-SMA+Snow-17(即薩克拉門托模型和Snow-17的組合模型)可以更好地完成日尺度上的降雨徑流模擬,能夠通過實測氣象數據進行流量預報;Le等[19]將LSTM模型用于湄公河受上游水電站和水庫泄洪影響流域的3 d預見期的洪水預報;Sahoo等[20]采用LSTM模型進行印度馬哈納迪河月尺度上的枯季流量預報,達到較高的預報精度;殷兆凱等[21]利用LSTM針對0~3 d的預見期分別建立了降雨徑流模型,說明在相同預見期下,LSTM模型擁有比新安江模型更好的預報能力。

本文以福建木蘭溪支流延壽溪渡里站以上流域為研究區,利用逐時實測降雨與流量數據建立BP神經網絡、LSTM神經網絡模型,進行未來1~24 h的逐時洪水預報,比較2種模型的預報結果,采用模塊化建模方法,討論了隱含層不同神經元數對模型訓練速度和預報精度的影響以及造成精度差異的可能原因,同時探討了如何簡單有效地避免模型訓練時陷入局部最優解,為將神經網絡模型應用于有水文站控制的小流域短期洪水預報工作提供參考。

1 研究方法

1.1 BP神經網絡模型

BP神經網絡通常包含3層結構,即輸入層、隱含層和輸出層,見圖1。BP神經網絡的構建過程主要分為2個階段:第一階段是信號的前向傳播,從輸入層經過隱含層,最后到達輸出層,計算出預測結果,并利用預測結果和真實值構成代價函數;第二階段是誤差的反向傳播,從輸出層到隱含層,最后到輸入層,利用梯度下降法,依次調節隱含層到輸出層的權重,輸入層到隱含層的權重。在反向傳播的過程中,根據誤差不斷調整優化各參數的值;不斷重復該過程,直到收斂。

圖1 BP神經網絡結構示意

1.2 LSTM神經網絡模型

注:×與⊙為矩陣元素積;+為相加。圖2 LSTM記憶單元細節結構

1.3 模擬結果的評價指標

從整體流量過程與單次洪水預報2個角度進行模型預報精度評價。以相對誤差RE、確定性系數DC作為誤差指標。其中:

(1)

(2)

確定性系數DC也稱Nash-Sutcliffe效率系數。其值在(-∞,1]之間,1表示預報結果完美擬合實測值。

根據GB/T 22482—2008《水文情報預報規范》對于單次洪水事件的預報結果,洪峰流量以實測洪峰流量的20%作為許可誤差;峰現時間以預報根據時間至實測洪峰出現時間之間距的30%作為許可誤差;徑流深以實測值的20%作為許可誤差,若該值大于20 mm,則取20 mm,若該值小于3 mm,則取3 mm。

當誤差小于規范所規定的誤差即可視為合格預報。合格預報次數與預報總次數之比的百分數為合格率,表示多次預報總體的精度水平。合格率按下式計算:

(3)

式中 QR——合格率;n——合格預報次數;m——預報總次數。

預報項目的精度按合格率或確定性系數的大小分為3個等級,見表1。

表1 洪水事件預報精度等級

2 研究流域與數據

2.1 研究區概況

實驗流域為福建莆田市木蘭溪支流延壽溪的渡里水文站以上集水區,流域內地形站點信息見圖3,流域面積68.57 km2,屬典型的亞熱帶季風氣候,年平均降雨在2 000 mm以上。流域地處山區,地形最大高差達842 m。流域內植被覆蓋良好,森林覆蓋率達69.3%。流域內無水庫調節,僅有少數塘壩,水體面積占0.02%。流域內有河道水文站1個(渡里),雨量站3個(渡里、下張隆和林兜)。渡里水文站以上集水區面積小,且不受水庫調節影響,能夠很好展現降雨徑流響應關系,且在福建省水文水資源勘測局的幫助下能夠獲得充足的降雨和徑流逐時數據,是進行小流域降雨徑流預報十分良好的研究區。

圖3 延壽溪渡里站以上集水區高程及站點分布

2.2 數據

本文選用來源于福建省水文水資源勘測局水雨情數據庫的渡里水文站逐時流量數據,渡里、下張隆和林兜的逐時降水數據,時間為2014年6月至2018年11月,采用泰森多邊形計算流域面平均雨量后作為降雨輸入。由于沒有2014年6月至2018年11月的降水預報數據,所以在進行未來第2~24 h逐時預報時用同期上一時刻的實測降雨數據代替降水預報數據。

利用MODIS 16 d尺度的植被指數(NDVI)數據(MODIS13A2),計算流域平均每16 d的多年平均植被指數,并采用三次樣條法插值成逐日數據,作為神經網絡反映流域內植被的季節變化的一項輸入。

3 模型構建與結果分析

3.1 神經網絡預報模型構建

神經網絡雖然已經在水文預報中得到了較為廣泛的應用,但是對建模環節中的一些問題仍然沒有系統性的解決方案,這些問題包括如何確定網絡結構(包括確定輸入節點,如何確定隱含層層數以及神經元數),如何確定訓練的迭代次數以及避免訓練過程中陷入局部最優解,以及如何考慮降雨徑流過程不同階段的動態特征差異等。

3.1.1模型結構的確定

a) 輸入節點。為便于比較,BP與LSTM模型將在輸入節點上將保持一致,采用同樣的輸入及預處理方法。2個神經網絡的降雨數據輸入節點通過計算當前流量與前期累積降雨的相關關系來確定。計算了1~6 h的逐小時降雨,6~12、12~18、18~24、24~30、30~36 h的6 h累計降雨,36~48、48~60、60~72 h的12 h累積降雨,結果見圖4。由圖4可見,1~6 h逐小時降雨,6~12、12~18、18~24 h的6 h累積降雨量與當前流量具有較高的相關性,相關系數都超過0.35,因此選取這9項累積降水數據作為降雨輸入。2個神經網絡的流量數據輸入節點數通過對流量數據進行了自相關(AR)和偏自相關(PAR)(圖5)分析確定。由圖5可見,流量具有很強的自相關特征(圖5a),即使30 h前的流量值與當前流量也具有較高的相關性,但是根據流量的偏自相關圖(圖5b),5 h之后的流量值與當前流量的相關系數小于0.2。因此,可以考慮將前1~5 h的流量值作為模型的輸入數據。但在進一步通過多次模型訓練嘗試后,確定流量數據滿足要求的最佳滯時為2,即將前1 h和前2 h流量值作為模型的輸入。除了上述的降水及流量數據輸入,考慮到植被截留過程會減少到達地面的有效降雨量,而流域內植被具有季節性變化,植被的年內生長變化對徑流量有著較大影響。所以將NDVI數據也作為神經網絡的一項輸入。

圖4 流量與前期不同時段累積降雨量的相關關系

a) 自相關

b) 偏自相關圖5 逐時流量數據的自相關和偏自相關分析

b) 輸入數據預處理。由于流量數據的數值尺度和單位與降雨量數據存在差異,以及激勵函數本身的特性,如sigmoid函數對于小于-5、大于5的輸入數據響應值不敏感,趨于0和1。所以在構建神經網絡模型時,必須對流量和降雨數據進行歸一化處理,使每個輸入、輸出節點的數據均值近似為0,方差為1。同時,為了縮小數據中極大值造成的影響,放大極小值的影響,在進行歸一化處理前先對數據取對數。為了避免數據中的0值影響對數計算,取0.01代替數據中的0值。經過歸一化處理后,數據的絕對值將變成一種相對關系,來消除量綱對模型運算的干擾[21],同時可加快梯度下降的收斂速度。

c) 隱含層層數及隱含節點數。神經網絡的建立過程中,隱含層層數和神經元數會對結構和訓練產生多方面的影響。增加神經網絡隱含層層數可以提高學習精度,Hornik等[22]證明一個擁有足夠多的神經元的單隱含層神經網絡就能夠完成從輸入數據到輸出數據之間任何可測量的函數關系,并達到期望的精度。但是同時讓模型結構變得更加復雜,也讓訓練時間大大增長。Villiers和Barnard[23]的研究表明,由2個隱含層組成的神經網絡的魯棒性較差,收斂精度也較低。所以在設計神經網絡時應優先考慮有1個隱含層的3層網絡。因此本文中BP和LSTM模型的隱含層層數都設定為1個。

對于隱含層內神經元數量,本文采用“試錯法”,通過比較不同神經元數量對驗證數據預報結果的影響,來確定合適的神經元數量。圖6為在不同神經元數下,以確定性系數DC(Nash效率系數)作為評價指標,BP模型與LSTM模型對未來1~24 h流量滾動預報精度的變化。

從圖6a、6b可以看出,BP模型的精度在n=10(圖6a黑線)時達到最高,Nash效率系數為0.987~0.711,當神經元數量繼續增加時,Nash效率系數會隨著神經元數量的增加而降低,同時訓練時間大大增長。所以BP模型最終神經元數量確定為10個。LSTM模型的精度隨著神經元數量的增加而提高,當n=64(圖6b黑線)時,未來1~24 h的滾動預報Nash效率系數最高,為0.969~0.747,當n>64后,Nash效率系數則出現了下降,同時訓練時間大大增長。所以LSTM模型最終神經元數量確定為64個。

從上述對比可見,BP神經網絡的最優隱含層節點數明顯小于LSTM的最優節點數據,其原因在于:不論是BP網絡還是LSTM網絡,過于復雜的神經網絡都可能帶來模型的過擬合問題,但這個過擬合問題對于BP網絡非常突出,而對于LSTM網絡,其本身的記憶與遺忘功能使得權重在更新時,不僅充分利用了當前時刻的數據特征,還利用了門的結構來記住還是遺忘前一時刻的數據特征,也就使得并非所有數據都對權重的更新產生影響,能夠一定程度上降低過擬合產生的可能,因而LSTM網絡的最優隱含層神經元數大大多于BP網絡。

a) BP模型

b) LSTM模型圖6 不同模型的神經元數對驗證數據預報精度的影響

3.1.2模型訓練

使用2014年6月至2016年6月的逐時流量數據、逐時降雨數據和NDVI數據,共3類分12個節點輸入模型,對模型進行訓練。

考慮到流域的降雨-徑流關系在不同的降雨量級以及不同的流量漲落階段可能具有一定的差異性,建立一個單一的全局模型無法讓模型很好地學習到不同階段流量漲落變化和降雨量的關系,因此本文在神經網絡建模中使用模塊化方法[24],首先按照閾值將降雨和流量數據劃分若干組,然后對每組數據分別訓練不同的模型。本文將降雨流量數據分成4組進行訓練:第一組為降雨大于1 mm且流量大于5 m3/s;第二組為降雨小于1 mm,流量大于5 m3/s;第三組為降雨大于1 mm,流量小于5 m3/s;第四組為降雨小于1 mm且流量小于5 m3/s。

模型訓練前,還需要確定包括模型的損失函數、學習率、批處理量參數(batch-size)、訓練集中的交叉檢驗數據比例(validation-split)和訓練最大迭代次數等。

本文中2種模型的損失函數均選擇均方誤差(MSE),計算公式如下:

(4)

學習率決定了模型權重更新的幅度,如果學習率很低,訓練會變得可靠,但是需要花費更多的優化時間;如果學習率很高,輸出誤差對權重的影響就越大,權重更新就越快,但可能無法收斂。但是最佳學習率的確定仍舊缺少最優的方法。本文中BP和LSTM模型的學習率依賴經驗分別定為0.20、0.01。

批處理量參數(batch-size)指模型訓練中梯度下降時每個批次所含的樣本數。該參數的大小影響模型的優化程度和速度。如果不設置該參數,意味著模型在計算時,是一次把所有的數據輸入模型中,然后進行梯度下降算法,由于在計算梯度時使用了所有數據,計算速度慢,需要更多的迭代次數來達到最優。而本文由于使用多年逐時數據進行建模,數據量較大,一次性把所有數據輸入模型,內存占用率高,會極大加重模型計算負擔,影響訓練速度和模型精度。所以經過試錯,2個模型的批處理量參數(batch-size)均設置為240。

訓練集中的交叉檢驗數據比例(validation-split)是指從數據集中切分出一部分作為驗證集,驗證集不參與訓練,并且在每次迭代時在驗證集中評估模型的性能,如計算損失函數。本文中2個模型的該參數設置為0.2。

為了確定2種模型的單次訓練的最大迭代次數,在上述參數設定的條件下,分別將BP和LSTM模型訓練一次,繪制2種模型的損失函數變化,見圖7。BP模型在迭代500次后損失函數的數值仍然高于LSTM模型100次迭代后損失函數的數值。所以為了避免訓練時間浪費的同時達到最佳訓練效果,確定BP模型最大迭代次數1 000次,LSTM模型迭代次數100次。

a) BP模型

b) LSTM模型圖7 不同模型損失函數隨迭代次數的變化

神經網絡的訓練過程對初始權重較為敏感。由于初始權重的不同,每次訓練都可能會得到一組不同的參數。為了獲得最優參數,一種方法是采用多組獨立初始權重值分別進行訓練,然后從多組訓練結果中選擇神經網絡模型[25],但這種最優網絡僅是對訓練數據最優,對未來未知的真實驗證數據未必最優。也就是說訓練中最佳的神經網絡模型對于未來的數據并不一定能達到最優的預測結果。Wang等[26]提出了一種簡單有效的解決局部最優的方法,即將一個模型訓練若干次(如20次),得到20組模型參數,選擇其中較好的10組,形成一組模型集合,并將這個模型集合的結果取平均值作為最終的輸出。本文采用了這種方法。

3.2 模型驗證

利用訓練好的模型做下一個小時(T+1)的流量預報,再以預報的T+1 h流量作為當前流量輸入,預報T+2 h流量,以此循環滾動,從而完成未來T+24 h的循環滾動預報。本文選擇2016年7月至2018年11月的預報結果,驗證模型預報精度。

a) 對整體流量過程的預報精度評價。BP與LSTM模型驗證期的24 h滾動預報結果的Nash效率系數見圖8。2個模型不同預見期的驗證期預報結果散點見圖9。由圖8可以看出,2種模型在第1 h的預報,Nash效率系數大于0.96(BP模型為0.975,LSTM模型0.968);隨著預見期的延長逐漸下降,但LSTM模型整體預報效果優于BP模型,并且預報精度的衰減速度大大慢于BP模型。LSTM模型滾動預報至第24 h的Nash效率系數仍達到0.74,達到乙級預報精度,而BP模型則降到0.51,符合丙級預報精度。由此可見,在整體流量過程預報值中,綜合而言,LSTM模型精度顯著優于BP模型。

圖8 BP、LSTM模型滾動預報結果Nash效率系數

a) LSTM模型1 h預見期

b) LSTM模型6 h預見期圖9 LSTM、BP模型不同預見期預報結果散點

c) LSTM模型24 h預見期

d) BP模型1 h預見期

e) BP模型6 h預見期

f) BP模型24 h預見期續圖9 LSTM、BP模型不同預見期預報結果散點

b) 對洪水事件的預報精度評價。挑選2016年7月后的15場洪水作為驗證集,將洪水事件期間的24 h滾動預報值分別與實測數據相比較,對2個模型對洪水事件中的預報精度進行了進一步驗證。根據GB/T 22482—2008《水文情報預報規范》,將15場洪水24 h內不同預見期內洪峰流量、峰現時間和徑流深3個評價指標的合格率繪制成折線,見圖10。其中,圖10a為洪峰流量誤差合格率,圖10b為峰現時間誤差合格率,圖10c為徑流深誤差合格率,圖10d為綜合3個誤差指標后模型的最終合格率,即當一場洪水3個指標都合格時,該場次洪水的預報即為一場合格洪水。LSTM模型3個指標的合格率以及綜合合格率都優于BP模型。其中LSTM模型洪水預報在預見期為24 h的預報合格率為60%,達到丙級預報精度,而BP模型為20%,未達到丙級項目精度等級。

a) 洪峰流量誤差合格率

c) 徑流深誤差合格率

d) 模型最終合格率續圖10 評價指標合格率折線

進一步以20160927和20180828 2場洪水為例,比較了BP模型與LSTM模型對洪水過程預報的能力。20160927場次洪水過程線見圖11,圖11a、11b、11c、11d預見期分別為1、6、12、24 h。在24 h的預見期內,LSTM和BP模型都可以較好反映洪水漲落過程。在預見期為1 h時,2個模型預報能力相當。隨著預見期的延長,BP模型出現了系統性地偏大,高估了整個洪水過程的洪量。20180828場次洪水過程線見圖12,圖12a、12b、12c、12d分別代表預見期為1、6、12、24 h。20180828場次洪水擁有2次漲水過程。在24 h的預見期內,2種模型都能夠反映出2次漲水退水過程。但是BP模型卻在2個洪峰處出現了嚴重的低估現象。在預見期6~24 h的預報結果中BP模型在第一個洪峰處的預報值高于了第二個洪峰的預報值,與實測數據不符,洪峰預報能力退化明顯。通過對2場洪水事件預報能力的對比,LSTM模型在預測洪峰流量以及洪水過程上擁有更好的性能,能夠很好地預報出洪水的漲退水過程以及洪量。

a) 預見期1 h

b) 預見期6 h

c) 預見期12 h

d) 預見期24 h圖11 預見期不同時20160927場次洪水預報效果比較

a) 預見期1 h圖12 預見期不同時20180828場次洪水預報效果比較

b) 預見期6 h

c) 預見期12 h

d) 預見期24 h續圖12 預見期不同時20180828場次洪水預報效果比較

4 討論

BP神經網絡與LSTM神經網絡模型都屬于黑箱模型,從本文應用結果來看,LSTM模型比BP神經網絡具有更好的逐時流量預報精度,其主要原因在于流域降雨徑流過程中的土壤水蓄量、地下水蓄量等狀態變量對河流流量具有巨大的影響,也就是說狀態信息對于預報精度的可靠性具有很大作用。LSTM神經網絡利用由遺忘門、輸入門和輸出門組成的記憶單元,來控制丟棄或增加信息,使得信息可以有選擇地通過,從而實現對時間序列過程狀態特征的記憶功能,而BP神經網絡缺乏這種對狀態信息的處理功能。

由于LSTM的記憶與遺忘功能使得累積誤差減少,這使得LSTM模型精度衰減速度遠遠小于BP模型,在循環滾動預報中,隨預見期的延長,精度的優勢不斷突顯。同時在本文應用中LSTM模型的隱含層神經元經過試錯確定為64個,BP模型設定為10個。在這種設置下,LSTM網絡在100次迭代次數下已經能夠使得損失函數達到最優,而BP需要近1 000次迭代才能帶到LSTM模型100次迭代的結果。如果將2個模型的隱含層神經元數都設定為64個,LSTM網絡的訓練速度將更明顯優于BP模型。

盡管LSTM網絡中的“門”結構使神經元之間能相互作用,可以使損失函數更好地收斂,使得LSTM趨近全局最優解的能力,LSTM模型仍然無法完全避免出現局部最優解的可能性。為此本文使用了Wang等[26]提出的解決方案,即訓練多次模型,計算結果后取平均的處理方式。該方法簡單易行,也確實在一定程度上避免了局部最優解對模型精度的影響。此外,還可以嘗試其他一些解決方案,比如使用變化的學習率或者蜂-蟻群算法等[27]方法,以期進一步提高預報模型的穩定性。

5 結論

本文使用BP和LSTM神經網絡分別構建了福建木蘭溪支流延壽溪小流域的降雨徑流預報模型,進行了預見期為1~24 h的逐時流量滾動預報,并對比了2個模型的預報精度,結果如下。

a) LSTM模型和BP模型在預見期為1 h時預報精度相當(BP模型為0.975,LSTM模型0.968),但隨著預見期的延長,LSTM預報精度的衰減速度大大慢于BP模型,BP模型24 h預見期的預報效率系數降至0.51,而LSTM模型為0.74,LSTM模型整體預報效果顯著優于BP模型。

b) 將降雨徑流過程按照一定閾值分開訓練的模塊化建模方法,使模型能更好地把握不同降雨階段以及洪枯階段的流量過程動態特征,從而有利于提高對流量過程的總體預報精度。

c) 對模型訓練多次然后從中取集合平均值的方法,可以在一定程度上消除模型訓練中參數優化過程陷入局部最優解的現象,提高模型預報精度。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 54pao国产成人免费视频| 日韩精品久久久久久久电影蜜臀| 国产成人无码Av在线播放无广告| 亚洲一级毛片免费观看| 日韩精品高清自在线| 久草中文网| 1769国产精品视频免费观看| 婷婷六月色| 国产乱人伦AV在线A| 国产精品白浆无码流出在线看| 久久婷婷人人澡人人爱91| 国产肉感大码AV无码| 亚洲国产中文在线二区三区免| 国产精品人成在线播放| 亚洲天堂久久| yy6080理论大片一级久久| 亚洲成在线观看| 国产精品主播| 人妻精品久久无码区| 奇米影视狠狠精品7777| 蜜芽国产尤物av尤物在线看| 伊人成色综合网| 国产男女XX00免费观看| 九色综合视频网| 色婷婷亚洲综合五月| 欧美一级色视频| 久久精品国产精品一区二区| 国产理论精品| 亚洲精品在线影院| 亚洲成人精品在线| 男人天堂伊人网| 91激情视频| 欧洲极品无码一区二区三区| 成人在线不卡视频| 成·人免费午夜无码视频在线观看| 国产不卡网| 性做久久久久久久免费看| 免费人成视网站在线不卡| 99青青青精品视频在线| 国产主播在线一区| 色婷婷狠狠干| 波多野结衣的av一区二区三区| 亚洲中文精品人人永久免费| 欧美日韩中文字幕二区三区| 亚洲无码37.| 国产成人久久777777| 黄色网址免费在线| 久热这里只有精品6| 日韩 欧美 国产 精品 综合| 国产久操视频| 国产精品v欧美| 青青草原国产精品啪啪视频| 色噜噜在线观看| 国产尤物视频在线| 波多野结衣无码AV在线| 一级一级特黄女人精品毛片| 91娇喘视频| 毛片免费试看| 美美女高清毛片视频免费观看| 先锋资源久久| 中国精品久久| 多人乱p欧美在线观看| 91一级片| 一级毛片免费不卡在线 | 午夜一级做a爰片久久毛片| 欧美专区在线观看| 国产高清在线精品一区二区三区| 国产精品视频999| 久久久久久久久18禁秘| 91精品日韩人妻无码久久| 国产免费怡红院视频| 毛片免费网址| 国产尹人香蕉综合在线电影| 国产精品美女自慰喷水| 中文字幕人妻无码系列第三区| 思思热精品在线8| 午夜性刺激在线观看免费| 精品精品国产高清A毛片| 在线观看精品国产入口| 日本久久网站| 免费在线a视频| 中文字幕在线日本|