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

基于不同ELM的西北旱區(qū)參考作物蒸散量模擬模型

2019-01-21 07:05:36張皓杰崔寧博胡笑濤龔道枝
中國農(nóng)村水利水電 2019年1期
關鍵詞:模型

徐 穎,張皓杰,崔寧博,2,馮 禹,胡笑濤,龔道枝

(1.四川大學 水力學與山區(qū)河流開發(fā)保護國家重點實驗室、水利水電學院,四川 成都 610065;2.南方丘區(qū)節(jié)水農(nóng)業(yè)研究四川省重點實驗室,四川 成都 610066;3.西北農(nóng)林科技大學 旱區(qū)農(nóng)業(yè)水土工程教育部重點實驗室,陜西 楊凌 712100;4.中國農(nóng)業(yè)科學院 農(nóng)業(yè)環(huán)境與可持續(xù)發(fā)展研究所作物、高效用水與抗災減損國家工程實驗室,北京 100081)

0 引 言

參考作物蒸散量(reference crop evapotranspiration,ET0)是計算作物需水量的關鍵因子,對灌溉系統(tǒng)設計和農(nóng)業(yè)水資源管理具有重要意義[1],其精準預報對作物需水量估算、水資源優(yōu)化調配及實現(xiàn)智慧灌溉等有重要價值。目前ET0計算模型有60多種,主要可歸為輻射法、溫度法、綜合法和水面蒸發(fā)法等[2]。1998年聯(lián)合國糧食與農(nóng)業(yè)組織(Food and Agriculture Organization of the United Nation,F(xiàn)AO)將Penman-Monteith(P-M)模型作為ET0計算的標準模型[3],P-M模型雖然是FAO推薦的最佳估算模型,但需要較完整的氣象資料,而很多地區(qū)由于標準氣象站點少且分布不均、設備不足等問題,常發(fā)生氣象資料缺失的情況,因此許多適用于較少氣象因子輸入的ET0計算模型被提出,如溫度法中的Hargreaves-Samani(HS)[4]、Blaney-Criddle[5]等模型,輻射法中的Priestley-Taylor(P-T)[6]、Makkink[7]等模型,經(jīng)驗公式法中的Irmark-Allen[8]等模型。

隨著機器學習(machine learning,MC)的發(fā)展,越來越多機器模型被應用于ET0模擬,如Kumar[9]、Gorka[10]等對BP神經(jīng)網(wǎng)絡(Back-Propagation,BP)模型在ET0模擬中的應用進行探討,證實BP神經(jīng)網(wǎng)絡在各自研究區(qū)域內模擬結果良好,但BP神經(jīng)網(wǎng)絡模型學習速率固定,收斂速度慢且存在局部極值的問題。Tabari[11]等采用支持向量機(support vector machine,SVM)和自適應模糊推理系統(tǒng)(adaptive neural fuzzy inference system,ANFIS)算法對ET0進行模擬,發(fā)現(xiàn)該模型精度高于Blaney-Criddle、Hargreaves、Priestley-Taylor等模型。崔遠來[12]等利用遺傳算法(genetic algorithm,GA)優(yōu)化神經(jīng)網(wǎng)絡(Artificial neural network,ANN)模型,克服了BP網(wǎng)絡輸入層、隱含層節(jié)點確定的盲目性,提高了進化神經(jīng)網(wǎng)絡(GA-ANN)的精度和適應性。張育斌[13]等基于典型的H-S模型、Mc Cloud模型和簡易神經(jīng)網(wǎng)絡算法提出五變量輸入、LM算法的三層網(wǎng)絡結構改進模型,發(fā)現(xiàn)改進模型在氣象資料欠缺、氣溫異常波動情況下的ET0模擬精度有明顯提升。

極限學習機自被提出已被廣泛用于多個領域。Abdullah[14]、馮禹[15,16]等分別在伊拉克地區(qū)和川中丘陵區(qū)采用極限學習機對ET0進行模擬,取得了良好的模擬效果,發(fā)現(xiàn)ELM模型運行速度快、精度高、泛化能力好。極限學習機在訓練過程中只需調整隱含層節(jié)點數(shù),便可獲得唯一最優(yōu)解,合理確定隱含層節(jié)點數(shù)對模型精度有重要意義。鄒偉東[17]等利用經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)方法確定隱含層節(jié)點數(shù),建立了日光溫室溫度濕度模擬模型,以改善因試湊法確定隱含層節(jié)點數(shù)引起的缺陷。但目前尚無人研究隱含層節(jié)點數(shù)對ELM模型模擬ET0精度影響。并且國內外學者在運用ELM模型模擬ET0時,均默認激活函數(shù)為“sig”,事實上還存在“sin”、“hardlim”、“radbas”等激活函數(shù),關于不同激活函數(shù)對于ELM模擬ET0精度的影響的研究尚待開展。

本文構建基于不同激活函數(shù)“sin”、“radbas”、“hardlim”的ELM模型來模擬ET0,分析其模擬結果提出適用于ELM模型的激活函數(shù),在精度較高的ELM模型中探究最優(yōu)隱含層節(jié)點數(shù),并同物理模型進行對比,驗證ELM模型精度,尋求適用于中國西北地區(qū)的ET0模擬模型,以期為中國西北旱區(qū)農(nóng)業(yè)水資源高效利用提供參考。

1 材料與方法

1.1 研究區(qū)概況

中國西北旱區(qū)位于31°35′~49°15′N,73°25′~110°55′E,包括陜西、新疆、甘肅、青海、寧夏、內蒙古西部等地,具有干旱缺水、荒漠廣布等特點。從旱區(qū)劃分標準來看,西北地區(qū)大部分地區(qū)降水在400 mm以下,多數(shù)地區(qū)在50~200 mm之間,屬于干旱、半干旱區(qū),而其蒸發(fā)量高達2 000~3 000 mm,更加劇了西北地區(qū)的干旱缺水。西北地區(qū)水資源總量為3 017.5 億m3,其中農(nóng)業(yè)用水量為2 492.5 億m3,占西北旱區(qū)水資源總量的82.6%,比全國農(nóng)業(yè)用水在水資源總量中占比高21.3%[18]。

本文選取西北地區(qū)若羌、綏德、烏魯木齊、西寧4個氣象站點(站點分布見圖1)1993-2016年逐日氣象數(shù)據(jù),包括日最高氣溫Tmax、日最低氣溫Tmin、日照時數(shù)n、距地面10 m高處的風速(計算時采用FAO風廓線關系[3]換算為2 m高度風速u2)和2 m高度相對濕度(relative humidity,RH)。氣象資料來源于國家氣象信息中心,經(jīng)過嚴格篩選核對,質量良好。研究取1993-2012年數(shù)據(jù)為訓練集,取2013-2016年數(shù)據(jù)為模擬集。

圖1 站點分布圖Fig.1 Distribution of meteorological station

1.2 研究方法

本文以1998年聯(lián)合國糧食與農(nóng)業(yè)組織(Food and Agriculture Organization of the United Nation,F(xiàn)AO)推薦的FAO-56 Penman-Monteith(P-M)模型計算結果作為ET0標準值,選取溫度法中的Hargreaves-Samani(H-S)、輻射法中的Makkink(MK)、Irmark-Allen(I-A)等3種在西北地區(qū)ET0計算精度較高的物理模型與P-M模型進行比較,各模型及計算公式如表1所示。

表1 參考作物蒸散量計算模型Tab.1 Calculation models of reference crop evapotranspiration

注:ET0為參考作物蒸散量,mm/d;Rn為凈輻射,MJ/(m2·d);G為土壤熱通量,MJ/(m2·d);Tmean為平均氣溫,℃;es為飽和水汽壓,kPa;ea為實際水汽壓,kPa;Δ為飽和水汽壓-溫度曲線斜率,kPa/℃;γ為濕度計常數(shù),kPa/℃;C0為轉換系數(shù),取0.000 939;Ra為大氣頂層輻射,MJ/(m2·d);Rs為太陽輻射量,MJ/(m2·d);α為經(jīng)驗系數(shù),取為1.26;n為日照時數(shù),h/d;RH為相對濕度,%。

1.3 模型構建

1.3.1 極限學習機原理

極限學習機(extreme learning machine,ELM)是一種基于單隱含層前饋神經(jīng)網(wǎng)絡(Single-hidden Layer Feedforward Neural Network,SLFN)的新算法,該算法由輸入層、隱含層和輸出層組成,ELM拓撲結構繪于圖2。

圖2 極限學習機(ELM)拓撲結構圖Fig.2 Topological structure of extreme learning machine

對于任意N個樣本(Xi,ti),假設輸入層神經(jīng)元個數(shù)為n,對應n個輸入變量(X1~Xn);隱含層有l(wèi)個神經(jīng)元(O1~Ol);輸出層神經(jīng)元個數(shù)為m,對應m個輸出變量(Y1~Ym);設輸入層與隱含層的連接權值ωl×n,隱含層與輸出層的連接權值βl×m;設隱含層神經(jīng)元的激活函數(shù)為g(ω,X,b) ,則ELM網(wǎng)絡的輸出表達式為:

(5)

式中:βj=(βj1,βj2,…,βjm)T為第j個隱含層神經(jīng)元與輸出變量的連接值向量;ωij=(ωj1,ωj2,…,ωjn)T為第i個輸入變量與第j個隱含層神經(jīng)元的連接值向量;bj為第j個隱含層神經(jīng)元的閾值。

由式(5)可得:

H·β=t

(6)

算法隨機產(chǎn)生輸入層與隱含層間的連接權值及隱含層神經(jīng)元的閾值且在訓練過程中保持不變,故只需設置隱含層神經(jīng)元個數(shù)l即可計算出隱含層與輸出層的連接權值,獲得唯一最優(yōu)解。因此隱含層與輸出層的連接權值β可以由式(6)的最小二乘解得:

(7)

式中:H+為H的MP廣義逆矩陣。

最小二乘解β求得后,訓練過程即完成。

根據(jù)ELM學習算法的步驟,可知其激活函數(shù)g(ω,X,b)與隱含層神經(jīng)元個數(shù)l對ELM模型的訓練和模擬有著顯著影響,故選擇合理的激活函數(shù)與合適的隱含層神經(jīng)元個數(shù)是至關重要的。

1.3.2 模型構建

僅設置ELM模型的隱含層節(jié)點數(shù)l與激活函數(shù)g(ω,X,b),根據(jù)初步模擬結果,設置隱含層節(jié)點數(shù)l為100,選取如下3個激活函數(shù)[19]:Sine函數(shù)、徑向基傳輸函數(shù)(radbas)、Hardlim函數(shù)。構建氣象因子輸入相同情況下不同激活函數(shù)的ELM模型。模型完整訓練過程及具體的程序見文獻[20]。

Sine函數(shù):

g(ωijX+bj)=sin(ωijX+bj)

(8)

徑向基傳輸函數(shù)(radbas):

(9)

Hardlim函數(shù):

(10)

1.4 模型驗證因子

選用6個評價因子[21]反映不同模型模擬精度,具體公式如下。

均方根誤差RMSE

(11)

決定系數(shù)R2:

(12)

納什系數(shù)NSE:

(13)

平均相對誤差MRE:

(14)

平均絕對誤差MAE:

(15)

整體評價指標[22]GPI:

(16)

2 結果分析

2.1 基于不同激活函數(shù)ELM模型的參考作物蒸散量日值模擬

分析基于“sin”、“radbas”、“hardlim”等3種激活函數(shù)的ELM模型在不同氣象因子組合輸入下的精度。ELM-sin、ELM-rad、ELM-hard的模擬結果分別如表2、表3和表4所示。

表2 不同氣象因子輸入下ELM-sin參考作物蒸散量日值模擬精度Tab.2 Reference crop evapotranspiration simulation accuracy of ELM-sin with different meteorological factors

注:** 表示1%水平上的極顯著相關,ELM-sin、ELM-rad、ELM-hard分別代表激活函數(shù)為sin、radbas、hardlim的ELM模擬模型,下同。

表3 不同氣象因子輸入下ELM-rad參考作物蒸散量日值模擬精度Tab.3 Reference crop evapotranspiration simulation accuracy of ELM-rad with different meteorological factors

表4 不同氣象因子輸入下ELM-hard參考作物蒸散量日值模擬精度Tab.4 Reference crop evapotranspiration simulation accuracy of ELM-hard with different meteorological factors

輸入5個氣象因子時,ELM-sin1模型與ELM-rad1模型的R2、NSE和RMSE分別為0.979 2和0.978 2、0.978 9和0.977 9、0.262 2 mm/d和0.269 4 mm/d,而ELM-hard1的R2、NSE和RMSE分別為0.935 2、0.933 6和0.484 8 mm/d。因此,ELM-hard1精度較低,其中ELM-sin1的可靠性最高。

輸入4個氣象因子時,ELM-sini、ELM-radj和ELM-hardk(i,j,k=2,3,4),GPI均排名前5,模型準確性較高。ELM-sini和ELM-radj(i,j=3,4)的R2和NSE均大于0.96,ELM-hardk(k=3,4)的R2和NSE均小于0.95,因此ELM-sin和ELM-rad模型較ELM-hard模型更好地反映了相同氣象因子組合與ET0的映射關系。并且,ELM-sin3和ELM-rad3的MRE均為8.96%,MAE分別為0.205 0 mm/d和0.203 8 mm/d,可見ELM-rad3較ELM-sin3模擬偏差更小;而ELM-sin4和ELM-rad4的R2和NSE相等,RMSE分別為0.327 0 mm/d和0.326 3 mm/d,ELM-sin4的RMSE更大,表明ELM-rad4準確性更高。同時,缺失u2的ELM-sin2、ELM-rad2和ELM-hard2模型較缺失RH、缺失n的ELM-sini、ELM-radj和ELM-hardk(i,j,k=3,4)模型精度更低,說明風速對西北旱區(qū)ET0的影響較大,這與汪彪[22]、劉憲峰[23]等得出的風速是影響西北地區(qū)的ET0計算主要因子之一的結論一致。

輸入3個氣象因子時,以u2、Tmax、Tmin作為輸入項的ELM-sin7、ELM-rad7和ELM-hard7的GPI均排名前4,三者的R2均大于0.90,其中ELM-hard7的R2和RMSE分別為0.941 0和0.453 6 mm/d,ELM-sin7和ELM-rad7的R2均大于0.96,RMSE分別為0.344 8和0.345 1 mm/d,說明ELM-sin7和ELM-rad7精度高于ELM-hard7,且ELM-sin7的模擬偏差更小。分析ELM-sini與ELM-radj(i,j=5,6,7)的模擬結果,u2對ET0模擬精度的影響最大,RH次之,n影響最小。Saeid[24]等在伊朗的研究表明,伊朗部分地區(qū)ET0對RH較為敏感,張青雯[25]等發(fā)現(xiàn)n是影響西南五省ET0的主導因素,而本研究中u2對ET0模擬精度的影響更為顯著,可見各氣象因子對不同地區(qū)ET0的影響和貢獻率有較大差異。

輸入2個氣象因子時,ELM-sini、ELM-radj和ELM-hardk(i,j,k=8,9)模型精度明顯下降。僅輸入溫度時,ELM-sin8的R2為0.797 1,而ELM-rad8和ELM-hard8的R2分別為0.782 8和0.787 7,說明僅輸入溫度時,ELM-sin8的精度最高,ELM-rad8次之,ELM-hard8精度最低。此外,結果表明基于溫度的ELM-sin8、ELM-rad8和ELM-hard8的精度較基于RH與u2的ELM-sin9、ELM-rad9和ELM-hard9精度更高,這與曹雯等[26]對西北旱區(qū)ET0敏感性分析得出的ET0對溫度輻射的敏感系數(shù)高于對相對濕度的敏感系數(shù)結論一致。

2.2 隱含層節(jié)點數(shù)對ELM模型精度的影響

僅輸入Tmax、Tmin、u2時,基于激活函數(shù)“sin”的ELM-sin7模型ET0模擬精度較高,可見該模型在氣象資料缺失時在西北地區(qū)適用性強,因此本文分析隱含層節(jié)點數(shù)目對ELM-sin7模型ET0模擬精度的影響。基于ELM-sin7模型,設置不同的隱含層節(jié)點數(shù),得出隱含層節(jié)點數(shù)與各模型精度評價指標的關系繪于圖3。

圖3 各模型驗證因子與隱含層節(jié)點數(shù)的關系Fig.3 Relationship between the verification factor of each model and the number of hidden layer nodes

由圖3可知,隨隱含層節(jié)點數(shù)的增加,RMSE呈先減小后趨于穩(wěn)定再增大的趨勢,MRE、MAE呈先減小后趨于穩(wěn)定再增大的趨勢,R2和NSE均呈先增大后趨于穩(wěn)定的趨勢。當隱含層節(jié)點數(shù)小于60時,模型模擬精度隨隱含層節(jié)點數(shù)增加明顯提高;當隱含層節(jié)點數(shù)達60~100時,RMSE、MAE、MRE、R2和NSE分別穩(wěn)定趨于0.345 2 mm/d、0.257 5 mm/d、10.8%、0.965 1和0.964 3;當隱含層節(jié)點數(shù)達到120后,RMSE與R2呈減小趨勢,模型精度開始降低。因此,隱含層節(jié)點數(shù)為60~100時ELM-sin7能取得較高的模擬精度。

增加隱含層節(jié)點數(shù)時,ELM模型能取得更高的精度,是由于增加隱含層節(jié)點數(shù)目可以實現(xiàn)將線性不可分的樣本映射到高層的線性可分空間,獲得更好的泛化能力。本研究僅選取了3種激活函數(shù)構建模型,基于輸入較少氣象因子仍能取得較高精度的ELM-sin7模型進行了隱含層節(jié)點數(shù)的討論,發(fā)現(xiàn)隨隱含層節(jié)點數(shù)增加,模型精度先有明顯提高后趨于穩(wěn)定,但當隱含層節(jié)點數(shù)過多時,出現(xiàn)精度下降的現(xiàn)象,說明隱含層節(jié)點數(shù)過多會導致過擬合問題,模型穩(wěn)定性遭到破壞。

2.3 ELM-sin在相同輸入因子下與物理模型精度比較

將基于溫度資料的H-S模型與ELM-sin模型進行比較,基于溫度和輻射資料的MK、I-A與ELM-rad模型進行比較,比較結果如表5所示。僅輸入Tmax、Tmin時,ELM-sin8的R2、NSE、RMSE和MAE分別為0.797 1和0.765 1、0.922 1和0.694 0 mm/d,H-S模型的R2、NSE、RMSE和MAE分別為0.762 7、0.594 2、1.199 4和0.915 2 mm/d,可見,ELM-sin8的精度明顯高于H-S模型。輸入Tmax、Tmin、n三個氣象因子時,ELM-rad5的R2、NSE、MAE、RMSE分別為0.822 5、0.785 8、0.664 1、0.882 3 mm/d,而I-A和MK的R2均小于0.77,NSE均小于0.63,RMSE均大于1.1 mm/d,MAE均大于0.9 mm/d,表明ELM-sin5的模擬精度明顯高于MK模型和I-A模型。因此,僅輸入溫度時推薦用ELM-sin8代替?zhèn)鹘y(tǒng)物理模型進行西北地區(qū)ET0模擬;僅輸入溫度和日照時數(shù)時,推薦用ELM-rad5模型。

表5 ELM-sin模型、ELM-rad模型與其他物理模型模擬精度比較Tab.5 Comparison of simulation precision between ELM-sin model、ELM-rad model and other physical models

2.4 ELM-sin7模型可移植性分析

基于激活函數(shù)“sin”的ELM-sin7模型(僅輸入Tmax、Tmin、u2)在較少參數(shù)輸入下ET0模擬精度較高,為檢驗ELM-sin7在西北旱區(qū)的普遍適用性,分別選取模擬站點P和訓練站點T,形成12組模擬集和訓練集樣本,構建ELM-sin7模型,模擬結果見表6。

表6 不同站點間ELM-sin7模型可移植性結果Tab.6 ELM-sin7 portability results among different stations

由表6可知,各站點組合下的ELM-sin7的R2均大于0.94,NSE均大于0.93,RMSE、MRE和MAE分別小于0.40 mm/d、14%和0.30 mm/d,模型精度較高,與站點組合前運行結果相比精度下降甚微,其中除以西寧作為訓練站點的3種站點組合外,其余的ELM-sin7的R2和NSE均大于0.96,MRE均小于12%,與原站點數(shù)據(jù)作為訓練集和測試集的模擬精度相當,且以綏德作為訓練站點的3種站點組合模擬精度較原有數(shù)據(jù)模擬精度明顯提高。因此,ELM-sin7在各站點之間的可移植性較好,在西北旱區(qū)的適用性強。

3 結 論

(1)基于“sin”、“radbas”、“hardlim”三種激活函數(shù)構建ELM模型,發(fā)現(xiàn)相同氣象因子輸入下,基于“sin”和“radbas”的ELM-sin和ELM-rad模型精度總是高于基于“hardlim”的ELM-hard模型,因此在模擬西北旱區(qū)ET0時不建議采用“hardlim”激活函數(shù)構建ELM模型。而氣象因子較少時,ELM-sin較ELM-rad能更精確地模擬ET0,其中以溫度和風速為輸入項的ELM-sin7模型的R2、NSE和RMSE分別為0.965 2、0.964 4和0.344 8 mm/d,因此,推薦基于溫度和風速的ELM-sin7作為西北旱區(qū)較少氣象因子輸入下的ET0模擬模型;僅以溫度作為輸入項時,ELM-sin8的RMSE、R2和NSE分別為0.922 1 mm/d、0.797 1、0.765 1,模型精度較ELM-rad8更高,因此,僅輸入溫度條件時,推薦基于溫度資料的ELM-sin8作為西北旱區(qū)的ET0模擬模型。

(2)對隱含層節(jié)點數(shù)進行研究表明,僅輸入Tmax、Tmin和u2的ELM-sin7模型的隱含層節(jié)點數(shù)小于60時,其精度隨隱含層節(jié)點數(shù)增加而提高;隱含層節(jié)點數(shù)為60~100時,ELM-sin7模型的RMSE、MAE、MRE、R2和NSE分別趨于0.345 2 mm/d、0.257 5 mm/d、10.8%、0.965 1和0.964 3;隱含層節(jié)點數(shù)大于120時,模型精度呈下降趨勢。因此,在僅輸入溫度和風速的情況下,利用ELM-sin7模型對西北旱區(qū)ET0進行模擬時,建議設置模型隱含層節(jié)點數(shù)為60~100。

(3)通過將ELM模型與H-S、MK和I-A模型進行對比發(fā)現(xiàn),僅輸入溫度時,ELM-sin8模型精度高于H-S模型;僅輸入溫度和日照時數(shù)時,ELM-rad5模型優(yōu)于MK模型和I-A模型。因此僅基于溫度或溫度和輻射資料進行西北旱區(qū)ET0模擬時,推薦使用ELM-sin8和ELM-rad5模型。

(4)ELM-sin7模型可移植性分析表明,在各訓練站點與模擬站點的組合下,ELM-sin7模擬精度較高,R2和NSE分別大于0.94和0.93,RMSE均小于0.40 mm/d。可見,在西北旱區(qū)內,ELM-sin7模型的可移植性和泛化能力較好,因此,某站點氣象資料缺失時,可用區(qū)域內臨近站點的氣象資料進行該站點的ET0模擬。

本文僅選擇了西北旱區(qū)的4個站點進行模型構建,在組合氣象參數(shù)時,參考其他文獻選取了精度最高的9種組合,未列出所有可能組合。在后續(xù)研究中應以區(qū)域實測蒸散量為模擬目標值,建立基于ELM的區(qū)域蒸散量模擬模型,以為西北旱區(qū)的高效節(jié)水灌溉和區(qū)域水資源管理提供更精確的科學依據(jù)。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 久996视频精品免费观看| 中文字幕亚洲精品2页| 国产精品专区第一页在线观看| 国产永久免费视频m3u8| 午夜福利无码一区二区| 国产精品午夜福利麻豆| 一本大道视频精品人妻| 色婷婷综合激情视频免费看| 日本www色视频| 真人高潮娇喘嗯啊在线观看 | 国产精品综合久久久| 欧美日韩在线亚洲国产人| 成人日韩精品| 少妇高潮惨叫久久久久久| 亚洲午夜18| 成人亚洲国产| 99激情网| 国产在线一区视频| 爆操波多野结衣| 久久国产免费观看| 国产超碰在线观看| 色综合久久久久8天国| 久久综合九九亚洲一区| 国产成人亚洲综合a∨婷婷| 精品国产网| 尤物在线观看乱码| 欧美日本在线| 日韩精品一区二区三区中文无码| 91无码视频在线观看| 亚洲网综合| 亚洲人成成无码网WWW| 国产欧美精品一区aⅴ影院| 青青久视频| 国产精品v欧美| 午夜无码一区二区三区| 欧美区在线播放| 亚洲精品视频免费| 久久青青草原亚洲av无码| 九色视频线上播放| 国产福利免费在线观看| 成人福利在线观看| 成人午夜亚洲影视在线观看| 色吊丝av中文字幕| 欧美亚洲中文精品三区| 五月婷婷综合网| 日韩二区三区无| 亚洲精品免费网站| 国产日产欧美精品| 蜜桃视频一区| 99久久国产精品无码| 一级毛片在线播放| 一区二区三区国产| 国产呦精品一区二区三区下载| 日本免费a视频| 人人看人人鲁狠狠高清| 国产激情无码一区二区APP| 亚洲精品综合一二三区在线| 国产成在线观看免费视频| 97se综合| 71pao成人国产永久免费视频| 亚洲第一视频网| 国产成人亚洲日韩欧美电影| 女人av社区男人的天堂| 日韩中文精品亚洲第三区| 亚洲专区一区二区在线观看| 亚洲av无码成人专区| 国产成人乱无码视频| 欧美日韩亚洲国产主播第一区| 99爱在线| 久久人人爽人人爽人人片aV东京热 | 欧美伦理一区| 亚洲无码免费黄色网址| 国产精品亚洲片在线va| 青青操视频免费观看| 久久青草精品一区二区三区| 五月天久久婷婷| 欧美成人午夜视频| 综合久久五月天| 98超碰在线观看| 亚洲成aⅴ人在线观看| 久久一色本道亚洲| 国产资源免费观看|