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

基于光能利用率模型的松嫩平原玉米單產估算

2017-09-12 07:36:47史曉亮楊志勇王馨爽薛羽梅
水土保持研究 2017年5期
關鍵詞:產量模型研究

史曉亮, 楊志勇, 王馨爽, 薛羽梅, 劉 峰

(1.西安科技大學 測繪科學與技術學院, 西安 710054; 2.中國水利水電科學研究院水資源研究所, 北京 100038; 3.國家測繪地理信息局陜西基礎地理信息中心, 西安 710054)

基于光能利用率模型的松嫩平原玉米單產估算

史曉亮1, 楊志勇2, 王馨爽3, 薛羽梅1, 劉 峰1

(1.西安科技大學 測繪科學與技術學院, 西安 710054; 2.中國水利水電科學研究院水資源研究所, 北京 100038; 3.國家測繪地理信息局陜西基礎地理信息中心, 西安 710054)

區域作物產量估測對于農業用地管理、糧食安全和宏觀調控具有重要意義。以中國重要的玉米主產區之一松嫩平原為研究區,利用2001—2013年長時間序列遙感和氣象數據,基于光能利用率模型估算玉米生育期5—9月累積地上生物量,進而根據地面實測的玉米收獲指數校正,建立了研究區玉米單產估算模型。結果表明:該模型估算結果達到了大范圍估產的精度要求;2001—2013年松嫩平原年均玉米單產為12 397.44 kg/hm2,總體呈增加趨勢。該區玉米單產空間分布差異顯著,總體為西低東高。從玉米單產年際變化的空間分布來看,研究區91.7%的區域玉米單產均有所增加,尤其是在松嫩平原西部的白城、齊齊哈爾等地,增產趨勢更為明顯。松嫩平原南部地區玉米單產與降水間具有顯著的正相關性,但與氣溫呈負相關性,而研究區其他地區玉米單產與降水和氣溫相關性不明顯。

松嫩平原; 玉米單產; 遙感; 光能利用率

由于人口的不斷增加和土地資源的快速減少,糧食安全問題在世界范圍內備受關注,因此,及其準確地掌握一個國家或區域的糧食產量及其動態變化,對于政府部門應對糧食安全問題,維護社會可持續發展具有重要意義[1-2]。傳統的作物產量估測多是采用人工區域調查方法,但不可避免地存在成本高、速度慢、統計主觀性強等問題,很難滿足大范圍估產的要求[3]。隨著遙感技術的發展,憑借其宏觀、及時和動態等特點,可以克服傳統方法的局限,為實現快速、準確的農作物估產提供了技術手段[4]。美國每年投入8 000多萬美元用于估計全球農作物產量,在世界糧食貿易中獲益18億美元[5]。中國于20世紀80年代后期開始作物長勢遙感監測和估產方面的研究[6]。經過近幾十年的發展,遙感估產取得了長足的發展,并得到了廣泛應用[7]。

目前,作物單產的遙感估算模型主要包括產量—遙感光譜指數的簡單統計相關模型、潛在—脅迫產量模型、產量構成三要素模型、作物干物質質量—產量模型等四類[7]。前三種模型均是通過將遙感信息及其輔助因子直接與產量相聯系構建模型,Balaghi等利用歸一化植被指數(NDVI)、降雨和溫度數據建立了摩洛哥小麥經驗回歸模型[8]。鄧坤枚等以內蒙古春小麥產區陳巴爾虎旗為研究區,利用春小麥乳黃熟期的環境減災衛星的NDVI數據和單產構建估產模型[9]。Ren等以山東省濟寧市為例,采用MODIS NDVI數據對冬小麥產量進行了估算[10]。但是此類模型均是針對特定區域,外推的適用性不高,而且缺少明確的生物物理機制,難以反映作物的生長生育過程。作物干物質—產量模型則是基于遙感數據估算地面干物質量,進而再依據其與果實間的關系得到作物產量。由于遙感數據一般主要反映作物整體的長勢狀況,相對于前三種模型直接估算產量,利用其估算干物質量更為合理,因此該方法成為當前作物單產估算的研究熱點,具有更加廣泛的適用性。目前計算作物干物質量較常用的方法是基于能量平衡的參數模型。任建強等[1,4]在黃淮海典型縣市基于CASA參數模型估算得到冬小麥關鍵生育期累積作物生物量,進而采用收獲指數加以校正,構建了作物生物量與產量間的關系模型,并證明了該方法估算冬小麥單產是可行的。

松嫩平原是中國重要的農產品生產基地之一,是北方春玉米帶的主要分布區,在保證國家糧食安全中起著至關重要的作用。因此,及時準確監測該地區的玉米產量,直接關系著國家農業經濟的健康發展,但是目前相關研究較少。本文擬利用2001—2013年長時間序列的遙感和氣象數據,基于光能利用率模型估算松嫩平原玉米生長季累積地上生物量,通過收獲指數校正估算得到研究區玉米單產,并驗證模擬精度,進而分析2001—2013年松嫩平原區域玉米單產的時空演變特征。以期為大范圍長時間序列作物估產提供一種可行的技術手段,并為相關政府部門制定糧食政策提供科學參考。

1 研究區概況

松嫩平原位于中國東北地區中西部,由松花江、嫩江沖積而成,是東北三大平原之一,也是我國重要的商品糧生產基地。松嫩平原南以松遼分水嶺為界、與遼河平原相隔,北與小興安嶺山脈相連,西以大興安嶺東麓丘陵和臺地為界,東至長白山外緣山麓臺地,是一個四周高、中部低、由周邊向中部緩慢傾斜的半封閉、不對稱的沉積盆地[11](圖1)。該區地處中國濕潤季風區與內陸干旱區之間的過渡帶,屬半干旱半濕潤氣候,春季干旱多風,夏季短促,冬季寒冷漫長。年均氣溫0~5℃,全年平均降水量400~600 mm,呈自西向東遞減趨勢[12]。松嫩平原的西部形成了農牧交錯區,中部為典型的農業耕作區,東部為低山丘陵區。

圖1松嫩平原在東北地區位置及氣象站點分布

2 材料與方法

2.1 數據來源與處理

MODIS NDVI數據來自美國NASA的地球觀測系統的2001—2013年玉米生長季5—9月的MOD13A3月值NDVI數據,空間分辨率為1 km。通過MODIS網站提供的專業處理軟件MRT(MODIS Reprojection Tools)對下載獲取的MODIS NDVI數據完成鑲嵌、裁切、投影與格式轉換等處理,得到松嫩平原長時間序列的MODIS NDVI數據集。

氣象數據來源于中國氣象科學數據共享服務網(http:∥cdc.nmic.cn),包括研究區19個氣象站點2001—2013年的逐月降水量、平均氣溫、逐日日照時數(圖1)。日太陽總輻射量根據童成立等建立的基于站點位置和日照時數的模擬方法估算得到[13],并匯總得到各氣象站月太陽輻射量。采用IDW插值方法,獲取與NDVI數據具有相同投影和像元大小的月太陽輻射、降水和平均氣溫柵格數據集。

2009年9月中下旬進行了玉米野外采樣,采樣范圍分布在整個研究區,隨機選擇32個樣點,利用收割法測量包括玉米、葉、莖等地上生物量,玉米單株產量,以及種植的株距和行距。選擇其中10個樣點估算研究區最大光能利用率和收獲指數,其余樣點實測數據用來對估算的玉米單產進行精度驗證。

2.2 作物單產估算方法

本研究采用光能利用率模型和遙感數據估算作物地上生物量,然后通過收獲指數校正地上生物量,估算得到作物產量[14-15],其計算公式如下:

Yield=B×HI

(1)

式中:Yield為作物產量(g/m2);B為作物生長季的累積地上生物量(g/m2);HI為作物收獲指數。

作物地上累積生物量可以由植被吸收的光合有效輻射(APAR)和光能利用率(ε)兩個因子來確定。

B=∑APAR×ε

(2)

式中:APAR表示作物吸收的光合有效輻射(MJ/m2);ε表示實際光能利用率(g/MJ)。作物吸收的光合有效輻射(APAR)計算方法為:

APAR=SOL×FPAR×0.5

(3)

式中:SOL表示太陽總輻射量(MJ/m2),可由大氣上界太陽輻射量和日照百分比計算。常數0.5表示植被利用的太陽有效輻射占太陽總輻射的比例,FPAR為植被層對入射光合有效輻射的吸收比例,計算公式為[16]:

(4)

式中:SRmin取值為1.08;SRmax的大小與植被類型相關,可由文獻獲取[17]。SR由NDVI計算得到。

(5)

光能利用率ε表示植被將吸收的光合有效輻射轉化為有機碳的效率(g/MJ),可由下式計算。

ε=Tε1×Tε2×Wε×εmax

(6)

式中:Tε1和Tε2分別表示低溫和高溫對光能利用率的影響系數,其中,Tε1為低溫和高溫情況下植物內在生化作用對光合的限制,由下式計算。

Tε1=0.8+0.02×Topt-0.0005×(Topt)2

(7)

式中:Topt表示研究區玉米生長季內NDVI值達到最大時當月的平均氣溫。Tε2表示環境溫度由最適宜溫度向高溫和低溫變化時植被光能利用率逐漸變小的趨勢。

(8)

式中:Wε為水分脅迫影響系數,表示水分條件對光能利用率的影響。隨著環境中有效水分的增加,其取值范圍為0.5(極端干旱條件)~1(非常濕潤條件),計算式為

Wε=0.5+0.5×EET(x,t)/PET(x,t)

(9)

式中:EET(x,t)為實際蒸散量(mm),可根據區域實際蒸散模型求取;PET(x,t)為潛在蒸散量(mm),可根據Boucher提出的互補關系求取[18]。εmax表示在理想條件下植被的最大光能利用率,其取值因植被類型不同有較大差別。對于收獲指數HI,在作物成熟期將作物地上部分進行實地脫粒、晾曬和稱重,將單位面積玉米粒質量與地上生物量之比確定為樣點的收獲指數。

2.3 數據分析方法

應用一元線性回歸分析方法在像元尺度分析松嫩平原2001—2013年玉米單產的變化趨勢。計算公式為[19]:

(10)

式中:θslope為趨勢斜率,即研究時段玉米單產的變化趨勢;n為研究時段的年數;Yieldi為第i年的玉米單產。斜率為正表示研究區玉米單產在13 a間的變化趨勢是上升的,反之表示下降。

此外,利用基于像元的相關分析方法,對2001—2013年松嫩平原玉米單產與主要氣候因子(降水與氣溫)進行逐像元相關分析,使用相關系數的空間分布情況來反映玉米單產與主要氣候因子的相關性。相關系數計算公式如下[20]:

(11)

3 結果與分析

3.1 最大光能率和收獲指數估算

在基于光能利用率的作物估產方法中,最大光能率和作物收獲指數是直接影響估算精度的重要參數。受到研究區所處地理位置、氣候條件、土壤性質等多種因素的影響,不同區域上述兩個參數的取值會有所不同。本研究中選擇10個樣點的實測玉米地上生物量與基于遙感數據估算的玉米生長季累積吸收的光合有效輻射建立線性模型,估算得到研究區最大光能利用率為3.08 g/MJ。在此基礎上,根據實測玉米產量和估算的地上生物量建立線性模型,獲取該地區玉米的收獲指數為0.46(圖2)。

圖2松嫩平原最大光能利用率和收獲指數估算

3.2 玉米單產估算精度驗證

為了驗證研究區玉米單產估算的精度,利用其余22個采樣點測量的單株玉米粒重及種植的株距和行距等參數計算樣點玉米單產,并與估算的玉米單產進行對比發現,松嫩平原2009年玉米預測單產的平均絕對誤差為-484.08 kg/hm2,平均相對誤差為-5.31%,基本達到了大范圍估產精度的要求。

3.3 松嫩平原2001-2013年玉米單產的時空變化特征

3.3.1 研究區玉米單產平均值的年際變化 由圖3可見,松嫩平原2001—2013年玉米單產年際波動變化較大,但整體呈增加趨勢。研究區玉米單產由2001年的10 423.77 kg/hm2增加到2013年的13 922.7 kg/hm2,年增速為167.55 kg/hm2(p<0.05),13 a間該區玉米平均單產為12 397.44 kg/hm2。從不同階段來看,2001—2002年,松嫩平原玉米單產增加顯著,之后到2008年玉米單產雖有波動,但整體呈增加趨勢,2007年產量較低,2008—2010年研究區玉米單產有所減少,此后呈顯著增加趨勢。

圖32001-2013年松嫩平原年均玉米單產變化動態

3.3.2 研究區年均平均單產的空間分布特征 松嫩平原玉米單產分布受到氣候、土壤類型、品種、耕作措施等多種因素的綜合影響,空間差異顯著,總體上呈西低東高。從圖4可以看出,在吉林省白城、扶余和黑龍江省的大慶等地,由于分布有范圍較大的鹽堿地,主要土壤類型為風沙土,土地生產力不高,且降水相對較少,因此玉米單產較低,多在8 000 kg/ hm2以下。而在研究區東部的黑龍江省的黑河市、哈爾濱市以及吉林省的長春市等地,主要土壤類型為黑土和草甸黑鈣土,有機質含量高,土壤肥沃,且氣候適宜,因此玉米平均單產大于13 000 kg/hm2。

圖42001-2013年松嫩平原年均玉米單產的空間分布

3.3.3 研究區年均平均單產變化趨勢的空間分布 本研究進一步采用逐像元線性趨勢法對2001—2013年松嫩平原玉米單產變化趨勢進行分析。從圖5可以看出,13 a來整個研究區的玉米單產發生了不同程度的變化,其中松嫩平原91.7%的區域玉米單產有所增加,尤其是在松嫩平原西部的白城市、齊齊哈爾市南部等地,由于自然條件較差,玉米產量較低,而隨著近年來節水增糧等項目的實施,玉米產量增加趨勢明顯。另有8.3%的區域玉米產量有所減少,主要分布在哈爾濱、長春市北部等部分地區,但減少趨勢不明顯。

3.4 研究區玉米單產與氣候因子的相關分析

降水和氣溫等氣候因子的變化對作物的生長發育具有重要影響,因此本文利用逐像元的相關分析方法,分別對2001—2013年研究區玉米單產、年降水量和年平均氣溫的相關系數進行分析,從而在空間尺度上分析松嫩平原玉米產量與主要氣候因子的關系,如圖6A所示,研究區玉米單產與降水量的相關系數為0.37,呈顯著正相關(p<0.05)。而從像元尺度來看,除了黑龍江省齊齊哈爾市的部分地區之外,松嫩平原大部分地區玉米單產與降水量之間表現為顯著正相關(p<0.05),呈顯著正相關的區域面積占全區玉米種植面積的比例達到55.69%(p<0.05)。尤其是在吉林省的白城、扶余和黑龍江省的大慶市南部等地玉米單產與降水量的相關性較高。

從圖6B可以看出,整個松嫩平原玉米單產與平均氣溫的相關系數為-0.16(p<0.05),不存在明顯的相關性。而從像元尺度來看,研究區玉米單產與氣溫呈顯著負相關的地區占玉米總種植面積的34.87%(p<0.05),主要分布在吉林省的白城、扶余和長春市,黑龍江省的南部地區,以及黑河市等地。研究區玉米單產與氣溫呈顯著正相關的地區占總面積的7.57%,主要分布在黑龍江省的綏化市、齊齊哈爾市等地。可以看出,松嫩平原大部分地區玉米單產與氣溫之間呈負相關性。

圖52001-2013年松嫩平原玉米單產變化趨勢

圖6松嫩平原玉米單產與年降水量、年平均氣溫間的相關性空間分布

4 結 論

本文利用長時間序列遙感、氣象等多源數據,在基于光能利用率模型計算松嫩平原玉米生育期累積生物量的基礎上,構建了研究區玉米單產估算模型,并分析了時空變化特征及其與氣候因子的關系,取得以下結論:

(1) 本研究基于作物干物質量—產量模式估算松嫩平原玉米單產,結果表明該方法的估算精度可以應用于大范圍作物估產研究。

(2) 松嫩平原2001—2013年平均玉米單產為12 397.44 kg/hm2,總體呈增加趨勢。玉米單產空間分布差異顯著,總體表現為西低東高。而從玉米單產變化趨勢的空間分布來看,研究區91.7%的區域玉米單產有所增加,尤其是在松嫩平原西部的白城市、齊齊哈爾市南部等地。

(3) 松嫩平原的吉林省白城、扶余、長春以及黑龍江省南部地區等地玉米單產與降水量具有較強的正相關,而與氣溫呈顯著負相關關系。而研究區的其他地區降水和氣溫對玉米單產影響不明顯。

[1] 任建強,陳仲新,唐華俊,等.基于植物凈初級生產力模型的區域冬小麥估產研究[J].農業工程學報,2006,22(5):111-117.

[2] 趙文亮,賀振,賀俊平,等.基于MODIS-NDVI的河南省冬小麥產量遙感估測[J].地理研究,2012,31(12):2310-2320.

[3] 馮美臣,肖璐潔,楊武德,等.基于遙感數據和氣象數據的水旱地冬小麥產量估測[J].農業工程學報,2010,26(11):183-188.

[4] 任建強,劉杏認,陳仲新,等.基于作物生物量估計的區域冬小麥單產預測[J].應用生態學報,2009,20(4):872-878.

[5] 吳炳方.全國農情監測與估產的運行化遙感方法[J].地理學報,2000,55(1):25-35.

[6] 錢永蘭,侯英雨,延昊,等.基于遙感的國外作物長勢監測與產量趨勢估計[J].農業工程學報,2012,28(13):166-171.

[7] 徐新剛,吳炳方,蒙繼華,等.農作物單產遙感估算模型研究進展[J].農業工程學報,2008,24(2):290-298.

[8] Balaghi R, Tychon B, Eerens H, et al. Empirical regression model using NDVI, rainfall and temperature data for the early prediction of wheat grain yields in Morocco[J]. International Journal of Applied Earth Observation and Geoinformation, 2008,10(4):438-452.

[9] 鄧坤枚,孫九林,陳鵬飛,等.利用國產環境減災衛星遙感信息估測春小麥產量:以內蒙古陳巴爾虎旗地區為例[J].自然資源學報,2011,26(11):1942-1952.

[10] Ren J Q, Chen Z X, Zhou Q B, et al. Regional yield estimation for winter wheat with MODIS-NDVI data in Shan-dong, China[J]. International of Applied Earth Observation and Geoinformation, 2008,10(4):403-413.

[11] 曾麗紅,宋開山,張柏,等.基于SEBAL模型與MODIS產品的松嫩平原蒸散量研究[J].干旱區資源與環境,2011,25(1):140-147.

[12] 黃志剛,王小立,肖燁,等.氣候變化對松嫩平原水稻灌溉需水量的影響[J].應用生態學報,2015,26(1):260-268.

[13] 童成立,張文菊,湯陽,等.逐日太陽輻射的模擬計算[J].中國農業氣象,2005,26(3):165-169.

[14] Bastiaanssen W G M, Ali S. A new crop yield forecasting model based on satellite measurements applied across the Indus Basin, Pakistan[J]. Agriculture Ecosystems and Environment, 2003,94(3):321-340.

[15] Moriondo M, Maselli F, Bindi M. A simple model of regional wheat yield based on NDVI data[J]. European Journal of Agronomy.2007,26(3):266-274.

[16] 朱文泉,潘耀忠,龍中華,等.基于GIS和RS的區域陸地植被NPP估算:以中國內蒙古為例[J].遙感學報,2005,9(3):300-307.

[17] 朱文泉,潘耀忠,張錦水.中國陸地植被凈初級生產力遙感估算[J].植物生態學報,2007,31(3):413-424.

[18] 羅玲.基于遙感—機理模型的松嫩西部草地凈初級生產力(NPP)模擬[D].長春:中國科學院東北地理與農業生態研究所,2010.

[19] 穆少杰,李建龍,周偉,等.2001—2010年內蒙古植被凈初級生產力的時空格局及其與氣候的關系[J].生態學報,2013,33(12):3752-3764.

[20] 史曉亮,李穎,趙凱,等.諾敏河流域植被覆蓋時空演變及其與徑流的關系研究[J].干旱區資源與環境,2013,27(6):54-60.

MaizeYieldEstimationBasedonLightEfficiencyModelinSongnenPlain,NortheastChina

SHI Xiaoliang1, YANG Zhiyong2, WANG Xinshuang3, XUE Yumei1,LIU Feng1

(1.CollegeofGeomatics,Xi′anUniversityofScienceandTechnology,Xi′an710054,China;2.DepartmentofWaterResources,ChinaInstituteofWaterResourcesandHydropowerResearch,Beijing100038,China;3.ShannxiGeomaticsCenterofNationalAdministrationofSurveying,MappingandGeoinformation,Xi′an710054,China)

Regional estimation of crop yield is critical for many applications such as agricultural land management, food security and macro-control. In the paper, the accumulated aboveground dry biomass in maize growing season from May to September every year were estimated using the long time sequenced remote sensing and meteorological data from 2001 to 2013 based on light use efficiency. Finally, we built the estimation model of maize yield through converting the calculated ground dry biomass to yield by harvest index. The results showed that the estimated yields can meet the precision requirements for estimating crop yield on a great scale, and it was feasible to predict maize yield by crop biomass based on remote sensing. Average annual maize yield from 2001 to 2013 was 12 397.44 kg/hm2in Songnen Plain, and the maize yield had increased over recent decades. There was significant spatial difference of maize yield, and the yield of the eastern part of the study area was higher than that of the western part. The spatial distribution of change trend of maize yield in study area from 2001 to 2013 was analyzed, which showed that about 91.7% of the study area presented an increasing yield trend, especially in the west of study area including Baicheng, Qiqihar, etc., where significant increase in maize yield had been seen. The maize yield had a significant positive correlation with the precipitation, and it had an obvious negative with temperature in the southern part of Songnen Plain.However, the maize yield had no obvious relationship with precipitation and temperature in the other area of the plain.

Songnen Plain; maize yield; remote sensing; light use efficiency

2017-01-17

:2017-02-23

水利部公益性行業科研專項經費項目(201401001);國家自然科學基金資助項目(51409204);西安科技大學科研培育基金資助項目(201605,2014008)

史曉亮(1985—),男,陜西寶雞人,講師,研究方向:資源環境遙感。E-mail:s_xiaoliang@126.com

楊志勇(1979—),男,湖南常德人,教授級高級工程師,主要從事水文水資源、分布式水文模擬,氣候變化對水資源影響等基礎研究。E-mail:yangzy@iwhr.com

S127

:A

:1005-3409(2017)05-0385-06

猜你喜歡
產量模型研究
一半模型
FMS與YBT相關性的實證研究
2022年11月份我國鋅產量同比增長2.9% 鉛產量同比增長5.6%
今年前7個月北海道魚糜產量同比減少37%
當代水產(2021年10期)2021-12-05 16:31:48
遼代千人邑研究述論
重要模型『一線三等角』
海水稻產量測評平均產量逐年遞增
今日農業(2020年20期)2020-11-26 06:09:10
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
主站蜘蛛池模板: 国产精品伦视频观看免费| 国模沟沟一区二区三区 | 伊人久久久久久久久久| 久久99热这里只有精品免费看| 97视频在线精品国自产拍| 欧美爱爱网| 国产嫩草在线观看| 国产一区自拍视频| 婷婷久久综合九色综合88| 91国内在线观看| 国产精品成人第一区| 久久伊人操| 一区二区三区国产精品视频| 福利国产在线| 五月天综合婷婷| 亚洲综合一区国产精品| 日韩国产精品无码一区二区三区| 日韩av电影一区二区三区四区| 国产精品自在自线免费观看| 久久精品电影| 永久天堂网Av| 日韩黄色精品| 伊人国产无码高清视频| 啪啪啪亚洲无码| 无码在线激情片| 伊人福利视频| 国产在线精彩视频论坛| 亚洲男人天堂久久| 波多野结衣无码视频在线观看| 国产精品无码AV中文| 久久香蕉国产线看观看精品蕉| 蜜桃臀无码内射一区二区三区| 国产一区二区视频在线| www.狠狠| 色屁屁一区二区三区视频国产| 国产免费怡红院视频| 免费不卡在线观看av| 九九热这里只有国产精品| 国产高颜值露脸在线观看| 国产成人在线小视频| 久久黄色小视频| 国产精品思思热在线| 91九色最新地址| 九九热精品在线视频| 欧美在线精品一区二区三区| 亚洲国产精品无码AV| 99久久99这里只有免费的精品| 日韩性网站| 色婷婷久久| 精品无码一区二区三区电影| 日本妇乱子伦视频| 亚洲成在线观看 | 91视频日本| 色综合久久88| 婷婷午夜影院| 国产精品9| 国产麻豆aⅴ精品无码| 热这里只有精品国产热门精品| 欧美日韩免费| 狠狠色香婷婷久久亚洲精品| 亚洲综合欧美在线一区在线播放| 一本久道久久综合多人| 高清免费毛片| 久久久精品无码一区二区三区| 婷婷丁香色| 精品丝袜美腿国产一区| 亚洲全网成人资源在线观看| 99在线视频精品| 亚洲国语自产一区第二页| 欧美国产精品不卡在线观看| 操国产美女| 国产成人精品一区二区三区| 精品综合久久久久久97| 九九热精品视频在线| 精品乱码久久久久久久| 国产精品太粉嫩高中在线观看| 99er精品视频| 成人日韩视频| 影音先锋丝袜制服| 久久亚洲黄色视频| 国产第八页| 制服丝袜 91视频|