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

采用作物系數法和PM模型估算南京地區玉米田蒸發蒸騰量

2016-03-23 02:56:53劉春偉吳心語邱讓建
節水灌溉 2016年9期
關鍵詞:模型

劉春偉,吳心語,邱讓建

(南京信息工程大學應用氣象學院,江蘇省農業氣象重點實驗室,南京 210044)

0 引 言

中國玉米種植面積位居世界第二,僅次于美國[1]。南京地區玉米種植水分供應來源主要是降水,玉米田蒸發蒸騰量大小的準確估算可以為提高玉米田水分利用效率提供參考,還可計算農田水量平衡和能量平衡,優化水資源配置。蒸發蒸騰量(Evapotranspiration,ET)的估算方法主要有直接法和間接法,直接法包括Penman-Monteith方法、Shuttleworth-Wallace方法等[2-4],間接法主要是單作物系數法和雙作物系數法[5,6]。隨著ET估算的發展,各種蒸散模型在各地區的適用性研究為農田水分管理提供了有利支持。

研究表明直接法和間接法基礎理論假設不同,實際應用上需要參照當地具體情況決定[7]。陳鳳等以大型蒸滲儀測定數據確定了陜西楊陵地區的夏玉米的作物系數,測定的夏玉米總ET為350 mm[8]。梁文清的研究表明作物系數法計算的夏玉米多年平均總ET約為340 mm[9]。眾多研究表明采用間接法估算夏玉米[10,11]、冬小麥的ET均取得良好的結果[12]。也有不少學者采用直接法估算冬小麥[4]、葡萄[13]等作物的ET,認為基于PM法的多源模型適宜估算稀疏作物的ET。直接法和間接法計算ET的方法各有優勢和缺點,在特定地區選用不同的模型計算ET有待具體分析。

本研究在測定氣象數據的基礎上,分別采用單作物系數法(Kc)、雙作物系數法(Kcb)和不同冠層阻力模型計算的PM方法(PM1和PM2)估算南京地區夏玉米ET,并以生育旺期液流法和微型蒸滲儀的實測ET為參考,擬得到南京地區夏玉米蒸散量的最優估算方法,然后分析整個生育期玉米田ET的變化,最后討論影響南京地區夏玉米蒸散量的關鍵氣象因素,為南京地區夏玉米田間水分管理提供理論依據。

1 材料與方法

1.1 試驗區概況

試驗地點位于南京信息工程大學農業氣象綜合實驗站(118.8°E,32.0°N,海拔32 m),試驗區屬亞熱帶季風氣候區,年平均降水量1 106 mm,年平均氣溫為15.6 ℃,年極端氣溫最高為39.7 ℃,最低為-13.1 ℃,土壤人為水成黃黏土,玉米品種為江玉403。該地區自然氣候資源豐富,年降水較多,適合夏玉米生長發育。

夏玉米株行距設置為50 cm×50 cm,于2015年6月18號播種,9月23日收獲,生育期劃分如下:6月18日到7月2日為播種出苗期;7月3日到8月8日為拔節期;8月9-19日為抽穗開花期;8月20-31日為灌漿期;9月1-18日為蠟熟期;9月19-23日為成熟期。

1.2 田間測量指標

(1)氣象數據。氣象要素的測定由Envis氣象站完成,測定指標包括2 m高度處每小時平均風速,m/s,氣溫,℃,相對濕度,%,凈輻射,W/m2,土壤熱通量,W/m2,降雨量,mm/d等。

(2)植株高度和葉面積指數的測定。植株高度通過鋼尺每周測定一次,葉面積指數通過照片法每周測定一次[4]。

(3)液流法測定植株蒸騰。在試驗中,使用液流法和微型蒸滲儀測定9月8-22日的夏玉米的實際蒸發蒸騰量。液流法中的莖部熱平衡法是指在莖或枝條外面裹上一個加熱套,連續加熱樹皮、木材和樹液,通過安裝在周圍的溫度傳感器來感應莖表面的溫度,依據熱量平衡原理求出被液流帶走的熱量來計算莖稈內液體的流量[14,15]。

(4)微型蒸滲儀測定土壤蒸發量。土壤蒸發采用微型蒸滲儀測定,蒸滲儀直徑D為10 cm,高h為20 cm。在夏玉米試驗田中不同位置布置3個微型蒸滲儀,每天18∶00采用精度為0.1 g的電子天平對微型蒸滲儀稱重,每兩次測定重量差值為ΔW(g),土壤蒸發E(mm/d)采用下式計算:

(1)

式中:ρ為水的密度,g/cm3。

1.3 單作物系數估算蒸發蒸騰量

單作物系數法其計算公式為[6]:

ETc=KcET0

(2)

其中:

(3)

式中:ETc為實際蒸發蒸騰量,mm/d;Kc為作物系數;ET0為參考作物蒸發蒸騰量,mm/d;T為2 m處的日平均氣溫,℃;u2為冠層上方2 m處的風速,m/s;Δ為飽和水汽壓與溫度關系曲線的斜率,KPa/℃;Rn為太陽凈輻射,MJ/(m2·d);G為土壤熱通量,MJ/(m2·d);es為飽和水汽壓,kPa;ea為實際水汽壓,kPa;es-ea為飽和水汽壓差,kPa;γ為濕度計常數,KPa/℃。

夏玉米不同階段的作物系數為Kcini=0.3,Kcmid=1.2,Kcend=0.6,因為當地的作物系數值需要按濕潤頻率和氣候條件對推薦值進行調整[6]。所以根據南京地區氣候條件,對Kcmid和Kcend進行修正,若日最低相對濕度的平均值不等于45%或2 m高度的日平均風速不等于2 m/s時,用公式(4)修正[6]:

(4)

式中:h為該生育階段內作物的平均高度,m;u2為2 m高度處的日平均風速,m/s;Rmin為日最低相對濕度的平均值,%;修正后,Kcmid和Kcend的值分別為1.14、0.58。

1.4 雙作物系數法

雙作物系數法是把作物系數Kc分為基礎作物系數和土壤蒸發系數兩個部分[6],基本作物系數Kcb用來描述作物蒸騰;另一個為土壤水蒸發系數Ke,反映降水或灌溉后地表土濕潤致使土壤蒸發強度短期內增加而對ETc產生的影響。公式為:

ETc=(Kcb+Ke)ET0

(5)

式中:ETc為作物需水量;ET0為參考作物蒸發蒸騰量,用公式(3)計算;Kcb為基礎作物系數;Ke為土壤蒸發系數。

基礎作物系數Kcb不同生育階段值按南京地區實際氣象資料根據公式(4)進行修正[6],Kcbini=0.15,Kcbmid=1.09,Kcbend=0.46。

當土壤表面由于灌溉或降水較濕潤時,Ke的值達到最大,但作物系數(Kc=Kcb+Ke)絕不會超過其最大值Kcmax,因為這個值是由土壤表層的騰發能量所決定的;當土壤表面干燥時,蒸發的水分較少,其蒸發的強度也小,這時可用下式計算[16]:

Ke=Kr(Kcmax-Kcb)≤fewKcmax

(6)

式中:Ke為土壤蒸發系數;Kr為表層土壤蒸發或累積深度的無量綱蒸發減小系數;Kcmax為降水或灌溉后的最大值Kc;Kcb為基礎作物系數;few為植株間發生蒸發土壤占總土壤面的百分比。

當土壤表面沒有可用蒸發的水分時,Ke最小,甚至可能為零。Kr的確定,當土壤表面濕潤(降雨1~2 d)時,Kr=1;降雨后(3~5 d),Kr=0.7,降雨后(6~8 d),Kr=0.2,表層土壤用來蒸發的總水量被消耗盡時,Kr=0。Kcmax[17]的計算公式為:

Kcmax=max({1.2+[0.04(u2-2)-

(7)

1.5 Penman-Monteith方法(PM1)計算蒸發蒸騰量

Penman-Monteith模型是將作物植株下墊面看作統一的整體來估算,又稱為單源模型,主要適用于下墊面均一,種植相對密集的作物[4]。計算方法如下:

(8)

式中:λ為汽化潛熱,2.45 MJ/kg;ρ為常壓下的空氣平均密度,kg/m3;Cp為空氣的比熱,1.013×10-3MJ/(kg·℃);rs是整體表面阻力,s/m;ra為空氣動力阻力,s/m。

空氣動力學阻力ra決定熱量和水汽由蒸發面傳輸到冠層以上大氣的過程,計算公式為[13,18]:

(9)

式中:Zm為測風速的高度,m;Zh為測濕度的高度,m;Zm=Zh=2 m;Zom為控制動量傳輸的粗糙長度,m;Zom=0.123h;Zoh為控制熱通量和水汽傳輸的粗糙長度,m;Zoh=0.1Zom;d為作物零面位移高度,m;d=2/3h。

總冠層阻力rs是對通過蒸騰和蒸發面水汽的阻力,當植被沒有完全覆蓋土壤表面,這個阻力必須包括土壤表面蒸發的影響。下式rs的計算如下[19]:

(11)

式中:c1=0.85;c2=1.83;ρa為大氣密度,kg/m3。

1.6 Penman-Monteith法(PM2)計算蒸發蒸騰量

冠層阻力rs采用Stannard(1993年)提出的基于飽和差、葉面積指數和太陽輻射的函數估算的總表面阻力,公式為[19]:

(12)

式中:e1=0.06,e2=1.00,e3=1 000;LAImax為日平均LAI的最大值;Rnmax為日平均太陽凈輻射的最大值,其他參數同上給出。

1.7 評價指標

為了衡量估算值與實際測定蒸發蒸騰量值的關系,本研究采用平均絕對誤差、一致性指數和決定系數來評價估算值,平均絕對誤差、一致性指數和決定系數的表達式如下:

(15)

這些指標直接反應了模擬值與實測值的相關關系及其絕對值大小。一致性指數(index of agreement,d)主要反映模擬值與實測值的離散程度和相對偏差[20]。d越接近1,MAE越小,表示模型估算值與實測值離散程度越小。R2的值越接近1,表示模型估算值與實測值相關性越顯著。

2 結果與分析

2.1 氣象數據生育期內的變化

圖1為2015年整個生育期氣象因素的變化。由圖1(a)可知,整個生育期日平均風速都維持在一個穩定的范圍,風速較小,最大值為2.42 m/s,出現在2015年8月9日。凈輻射值在整個生育期都較大,晴天條件下,凈輻射呈增大的趨勢,最高值在12~15 MJ/(m2·d)范圍內,陰雨天氣,出現低峰值,相對濕度較大。但從圖1(b)中可以看出,整個生育期降雨較少,相對濕度一直較低。空氣溫度在整個夏玉米生育期的波動不大,在8月初出現最高溫度32.5 ℃,隨后逐漸減小,在9月13日減小到最低值19.6 ℃。從圖1(c)可以看出,夏玉米整個生育期內的總降水量較少,降水量最大值出現在6月27日,次大值出現在8月10日,土壤水分含量隨降水的增加而增加,在降水過后又恢復到一個相對穩定的范圍,在7月下旬,隨著天氣變熱,凈輻射出現最大值,土壤含水量相應地出現最小值,土壤蒸發變大,這與南京地區7、8月的高溫天氣有關。

圖1 氣象因素生育期內的變化

2.2 模型估算值與液流值的比較

將9月8-22日時段單作物系數法Kc、雙作物系數Kcb、PM模型(PM1和PM2)估算夏玉米蒸發蒸散量的值與液流法和微型蒸滲儀的實測值進行比較分析(圖2和表1)。結果表明,采取作物系數法和PM模型估算玉米蒸發蒸散量的估算值與實測值變化趨勢總體一致,估算值相對比較穩定。從圖2和表1中可看出,單作物系數法估算的ET決定系數R2為0.57,平均絕對誤差MAE為1.28 mm/h,一致性指數d為0.33,其估測值較符合實測值,離散程度小;圖2(b)為雙作物系數法估算的ET與實測值間的關系,雙作物系數模型中R2為0.42,MAE為1.28 mm/h,d為0.35;圖2(c)為PM1模型估算值ET與實測值的關系,R2為0.52,MAE為0.79 mm/h,d為0.48,二者的顯著關系強于其他模型,建議使用PM1模型估算南京地區夏玉米的蒸發蒸騰量;圖2(d)為PM2模型估算的ET與實測值的關系,其擬合性最差,顯著性不明顯,R2為0.33,MAE為1.67 mm/h,d為0.28;表1給出了2015年南京地區夏玉米4種模型估算值的平均值,PM1模型估算的日平均ET值最大,為3.16 mm,與實測值的日平均值相等;其次是PM2模型和雙作物系數法估算的日ET,分別為2.93和2.55 mm;單作物系數法估算的日平均ET值較小,與其他模型的估算值相差顯著,為2.12 mm。綜合表明,采用PM1模型估算南京地區夏玉米蒸發蒸騰量最貼切,計算誤差較小。通過R2和MAE的比較,本研究采用PM1模型估算的夏玉米ET。

圖2 液流法和蒸滲儀測定夏玉米蒸散量和4種方法估算ET的關系

模型模型平均值/(mm·d-1)實測平均值/(mm·d-1)AbR2MAEdKc2.12Kcb2.55PM13.16PM22.933.160.65-0.200.571.290.330.64-0.100.421.290.350.630.380.520.800.480.68-0.590.331.680.28

注:線性公式為y=ax+b,a為直線的斜率,b為直線的截距,R2為決定系數,MAE為平均絕對誤差,d為一致性指數。

2.3 全生育期ET的變化

由圖3可知,夏玉米蒸發蒸騰量在生育期中期較大,生育后期減小,最大值出現在抽穗期或灌漿期,最小值出現在播種出苗期;從種植夏玉米開始,南京進入初夏,隨著太陽輻射和大氣溫度的增加,夏玉米的參考作物蒸發蒸騰量一直增加,在抽穗期達到最大值,到灌漿后期,太陽輻射與大氣溫度逐漸減小從而使收獲期的夏玉米參考作物蒸發蒸騰量降低。夏玉米作物蒸發蒸騰量主要集中在拔節期和抽穗灌漿期,整個生育期的變化趨勢呈單峰型[20]。從表2中可看出,南京地區夏玉米日平均蒸發蒸騰量為3.16 mm/d;播種出苗期ET為40.6 mm/d,占全生育期的13%,日平均蒸散量為2.71 mm/d;拔節期ET占44%,日平均蒸散量為3.7 mm/d;抽穗開花期ET占9%,日平均蒸散量為2.57 mm/d;灌漿期ET占14%,日平均蒸散量為3.59 mm/d;蠟熟期ET占17%,日平均蒸散量為2.88 mm/d;成熟期ET占3%,日平均蒸散量為2.08 mm/d;從模型估測值的日平均ET可看出,它們與液流法和微型蒸滲儀實測的日平均值的誤差≤0.34,誤差在合理范圍內,說明估測值在一定的程度上符合實際蒸發蒸騰量,并支持理論依據。

圖3 PM1法估算的玉米田ET的生育期變化

2.4 ET與氣象因素的關系

表2 PM1法計算不同生育期的ET

氣象環境因素是影響夏玉米蒸發蒸騰量的大小的主要因素,南京地區位于長江中下游,夏季高溫炎熱,冬季寒冷干燥,2014年年平均氣溫15.9~17 ℃,夏季平均氣溫25.4~26.2 ℃[21]。太陽輻射是地球能量的來源,它影響著作物體內水的汽化,大氣和土壤剖面的增溫,隨著全球天氣變暖的趨勢,夏玉米的生長發育在一定程度上受到了影響,但風速的大小與凈輻射的變化究竟哪一個對夏玉米蒸發蒸散量的影響更大呢?如何有效地應對氣象條件變化對蒸發蒸騰量的影響呢?圖4為PM-1模型估算的ET與氣象因素的關系,通過分析它們之間的相關性,可以明確ET對氣象因素的響應,從而為南京地區夏玉米種植提供依據。

圖4 ET與氣象因素的關系(注:除風速外,顯著性水平p<0.001)

從圖4可看出,夏玉米的蒸發蒸騰量與風速、凈輻射、溫度、飽和水汽壓差的決定系數R2分別為0.02、0.97、0.34和0.65,其中風速與ET呈不明顯的負相關關系,與凈輻射、溫度和飽和水汽壓差存在顯著的正相關,其相關性由大到小排序為凈輻射>飽和水汽壓差>溫度。此外,蒸發蒸騰量還受土壤條件和作物的生物學特性影響,強小嫚等(2008年)就氣孔阻力與ET的關系進行深入研究,發現在土壤含水率變化梯度較小時,作物蒸發蒸騰量在一定的氣孔阻力范圍內隨著氣孔阻力的增大而減小,當土壤含水率變化梯度較大時,氣孔阻力與蒸發蒸騰量的關系不明顯,原因是灌水下限較低時,土壤含水率變化快,使得田間水分狀況變復雜,從而影響了作物蒸發蒸騰[21]。

3 討 論

本研究采用單作物系數法、雙作物系數法、PM1和PM2方法估算了ET,并與液流微型蒸滲儀法測定ET進行比較,結果表明采用PM1模型估算的ET最貼近實際蒸發蒸騰量。盧曉鵬等運用云南曲靖市陸良站1990-1992年的逐日氣象資料對夏玉米蒸發蒸騰量進行估算,結果表明單作物系數法和雙作物系數法估算的ET與實測值ET非常接近,相對偏差都小于10%,相比而言,雙作物系數法估算的ET更切合實際,這主要是因為雙作物系數法充分考慮了土壤蒸發和植株蒸騰在各階段的變化規律,同時還考慮了灌溉和降水對作物系數的影響[22]。冠層阻力公式選擇也會影響估算的蒸發蒸騰量。丁加麗等(2003-2004年)運用PM模型和改進的PM模型估算江西省鷹潭余江示范區的稻田蒸散量[2,3],改進后的PM模型能反映出不同土壤水分狀況下冠層阻力的變化規律,與實測值ET的誤差從18.57%降低至10.84%,模擬效果較好,且與氣象因子的敏感性分析響應程度為凈輻射>平均溫度>風速,這與本文PM1模型估算玉米蒸發蒸騰量的情況相類似,其絕對誤差最小,能夠很好地模擬夏玉米蒸發蒸騰,且ET對氣象因子的響應順序基本一致。

4 結 語

本研究在試驗數據基礎上對南京地區夏玉米ET進行估算,并分析了ET對氣象因素的響應。ET估算方法包括單作物系數法、雙作物系數法、PM1和PM2方法,并以液流法和微型蒸滲儀的實測ET為參照。結果表明,相比于單作物系數法、雙作物系數法和PM2法,PM1法估算的ET精度更高,與實測值更貼近,其R2為0.52,MAE為0.8 mm/h,d為0.48。采用PM1法估算的全生育期ET,發現夏玉米ET主要集中在拔節期和抽穗灌漿期,其整體變化呈單峰型。ET對氣象因素的響應分析表明凈輻射對ET的影響最為顯著,其次是飽和水汽壓差和溫度,最后是風速。通過PM1模型對南京地區夏玉米ET的模擬可以為南京地區夏玉米的農田水分管理提供科學依據,從而制定有效的田間水分管理計劃。

[1] 何奇瑾,周廣勝. 我國玉米種植區分布的氣候適宜性[J]. 科學通報,2012,(4):267-275.

[2] Ershadi A, McCabe M F, Evans J P, et al. Impact of model structure and parameterization on Penman-Monteith type evaporation models[J]. Journal of Hydrology, 2015,525,521-535.

[3] 劉春偉,邱讓建,王振昌,等. 基于液流量的蘋果樹蒸騰量模擬[J]. 農業機械學報,2016,(2):105-112.

[4] 劉春偉,曾勰婷,邱讓建. 用分時段修正雙源模型估算南京地區冬小麥生育期蒸散量[J]. 農業工程學報,2016,(S1):80-87.

[5] Allen R G, Pereira L S, Smith M, et al. FAO-56 dual crop coefficient method for estimating evaporation from soil and application extensions[J]. Journal Of Irrigation And Drainage Engineering-Asce, 2005,131(1):2-13.

[6] Allen, R.G., Pereira, L.S., Raes, D. and Smith, M., 1998. Crop evapotranspiration[R]. FAO irrigation and drainage paper No. 56, No.56.

[7] 許 迪,劉 鈺. 測定和估算田間作物騰發量方法研究綜述[J]. 灌溉排水, 1997,(4):56-61.

[8] 陳 鳳,蔡煥杰,王 健,等. 楊凌地區冬小麥和夏玉米蒸發蒸騰和作物系數的確定[J]. 農業工程學報,2006,(5):191-193.

[9] 梁文清. 冬小麥、夏玉米蒸發蒸騰及作物系數的研究[D].陜西楊凌:西北農林科技大學,2012.

[10] 盧曉鵬,段順瓊,馬顯瑩,等. 單雙作物系數法計算玉米需水量的對比研究[J]. 節水灌溉, 2012,(11):18-21.

[11] 趙麗雯,吉喜斌. 基于FAO-56雙作物系數法估算農田作物蒸騰和土壤蒸發研究----以西北干旱區黑河流域中游綠洲農田為例[J]. 中國農業科學, 2010,(19):4 016-4 026.

[12] 樊引琴,蔡煥杰. 單作物系數法和雙作物系數法計算作物需水量的比較研究[J]. 水利學報,2002,(3):50-54.

[13] Zhang B, Kang S, Li F, et al. Comparison of three evapotranspiration models to Bowen ratio-energy balance method for a vineyard in an arid desert region of Northwest China[J]. Agricultural and Forest Meteorology, 2008,148(10):1 629-1 640.

[14] Chabot R, Bouarfa S, Zimmer D, et al. Evaluation of the sap flow determined with a heat balance method to measure the transpiration of a sugarcane canopy[J]. Agricultural Water Management. 2005,75(1):10-24.

[15] Logsdon SD, Singer JW, Prueger JH, et al. Comparison of corn transpiration, eddy covariance, and soil water loss[J]. Soil Science Society of America Journal. 2014,78(4):1 214-23.

[16] 范曉慧. 不同水分脅迫下青貯玉米需水量及優化灌溉制度的分析研究[D].呼和浩特:內蒙古農業大學,2013.

[17] 閆浩芳. 內蒙古河套灌區不同作物騰發量及作物系數的研究[D].呼和浩特:內蒙古農業大學,2008.

[18] Katerji N,Rana G,Fahed S. Parameterizing canopy resistance using mechanistic and semi-empirical estimates of hourly evapotranspiration: critical evaluation for irrigated crops in the Mediterranean[J]. Hydrological Processes, 2011,25(1):117-129.

[19] Li S, Zhang L, Kang S, et al. Comparison of several surface resistance models for estimating crop evapotranspiration over the entire growing season in arid regions[J]. Agricultural and Forest Meteorology, 2015,208:1-5.

[20] Moriasi D, Arnold J, Van Liew M, et al. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations[J]. Transaction of ASABE, 2007,50(3):885-900.

[21] 強小嫚.ET0計算公式適用性評價及作物生理指標與蒸發蒸騰量關系的研究[D]. 陜西楊凌:西北農林科技大學,2008.

[22] 盧曉鵬,段順瓊,馬顯瑩,等. 單雙作物系數法計算玉米需水量的對比研究[J]. 節水灌溉,2012,(11):18-21.

[23] 丁加麗,彭世彰,徐俊增,等. 基于Penman-Monteith方程的節水灌溉稻田蒸散量模型[J]. 農業工程學報, 2010,26(4):31-35.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 拍国产真实乱人偷精品| 国产精品手机视频一区二区| 久久夜色精品| 又爽又大又黄a级毛片在线视频 | 成人av专区精品无码国产| 欧洲日本亚洲中文字幕| 婷婷色在线视频| 素人激情视频福利| 精品视频免费在线| 呦女亚洲一区精品| 青青草国产一区二区三区| 欧美一区精品| 成人小视频在线观看免费| AV网站中文| 国产精品亚欧美一区二区| www.91中文字幕| 在线播放精品一区二区啪视频| 亚洲第一视频免费在线| 沈阳少妇高潮在线| 无码日韩精品91超碰| 日韩高清一区 | 亚洲中文字幕97久久精品少妇| 大香伊人久久| 国产精品嫩草影院av | 成人日韩欧美| 一级毛片免费高清视频| 综合亚洲色图| 久久免费视频6| 国产真实自在自线免费精品| 超碰免费91| 一级不卡毛片| AV片亚洲国产男人的天堂| 国产成人凹凸视频在线| 亚洲第一精品福利| 草逼视频国产| 亚洲综合色婷婷中文字幕| 日本高清免费一本在线观看 | 91免费观看视频| 青青草原国产| 精品欧美一区二区三区久久久| 国产精品视频观看裸模| 国产欧美日韩视频怡春院| 色综合中文| 中字无码av在线电影| 国产精品亚欧美一区二区三区 | 手机看片1024久久精品你懂的| 国产原创演绎剧情有字幕的| 毛片最新网址| 又爽又大又光又色的午夜视频| 伊伊人成亚洲综合人网7777| 国产一区二区三区精品久久呦| 麻豆精品久久久久久久99蜜桃| 99视频在线免费看| 在线播放91| 亚洲资源站av无码网址| 免费福利视频网站| 亚洲AV成人一区国产精品| 亚洲毛片在线看| 成人在线亚洲| 亚洲最猛黑人xxxx黑人猛交| 91黄色在线观看| 在线观看国产精品第一区免费 | 国产sm重味一区二区三区| 久久香蕉国产线看观看精品蕉| 国产簧片免费在线播放| 热久久这里是精品6免费观看| 日本尹人综合香蕉在线观看 | 日韩精品少妇无码受不了| 亚洲中文在线看视频一区| 亚洲色成人www在线观看| 91九色最新地址| 国产成人亚洲欧美激情| 精品91自产拍在线| 欧美在线观看不卡| 精品一區二區久久久久久久網站 | 国产在线拍偷自揄观看视频网站| 国产99在线| 热思思久久免费视频| 亚洲三级影院| 在线国产毛片| 国产 日韩 欧美 第二页| 无码人中文字幕|