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

內蒙古羊草草原物候及其對氣候變化的響應

2016-06-05 14:57:58芹,趙勝*,鄭
地理與地理信息科學 2016年6期

范 德 芹,趙 學 勝*,鄭 周 濤

(1.中國礦業大學(北京) 地球科學與測繪工程學院,北京 100083;2.北京師范大學資源學院,北京 100875)

內蒙古羊草草原物候及其對氣候變化的響應

范 德 芹1,趙 學 勝1*,鄭 周 濤2

(1.中國礦業大學(北京) 地球科學與測繪工程學院,北京 100083;2.北京師范大學資源學院,北京 100875)

為了有效分析內蒙古羊草草原的物候及其對氣候變化的響應,基于NOAA歸一化差值植被指數(NDVI)時序數據,利用S-G濾波擬合方法對原始NDVI時序數據進行了重建,然后采用滑動平均法對重建后的NDVI數據進行了物候期識別,利用地面物候觀測數據對識別結果的驗證結果表明,“S-G濾波+滑動平均法”識別的物候期與地面物候觀測結果具有較好的一致性。在此基礎上,獲取了各氣象站點的物候期及其生長季長度,結合氣候數據進行物候期和生長季長度對氣候變化的響應分析顯示,返青期呈提前趨勢,主要與當年春季溫度升高和前一年10月到當年4月的累積降水量增加相關;黃枯期呈提前趨勢,主要與當年8月份溫度升高相關;生長季長度呈延長趨勢,主要是由返青期提前的幅度較黃枯期提前幅度更大所致。

羊草;草原;NDVI;物候;遙感;氣候

0 引言

內蒙古草原植被是中國北方溫帶草原的主體,具有極其典型的代表性[1],羊草(Leymuschinensis)是其主要植被類型之一[2],對氣候變化的響應敏感。因此,研究羊草草原物候及其對氣候變化的響應,不僅有助于增進內蒙古草原植被對氣候變化響應的理解,還有助于研究全球植被生產力[3]、群落組成及結構、土壤-植被-大氣系統水熱碳交換的變化[4,5],對提高氣候-植被之間物質與能量交換的模擬精度、準確評估植被生產力與全球碳收支具有重要意義[6,7]。

目前,關于內蒙古草原植被物候及其對氣候變化響應的監測,主要采用地面觀測和遙感監測兩大類方法。地面觀測方法主要是根據內蒙古各地區多年的農氣觀測站資料和氣候數據,統計分析內蒙古典型草原植物物候對氣候變化的響應[8-11]。相關研究表明,近30年來內蒙古典型草原的物候期呈提前趨勢[8],其返青期的提前與春季氣溫升高呈正相關,與日照時數呈負相關,與降水的相關性因地域而異;黃枯期提前與前1~2個月溫度呈負相關,與日照時數呈正相關[10,12]。地面觀測法主要是通過人工觀測記錄植物的物候期[13],在物種水平上進行物候對氣候變化的響應研究[14],很難從群落甚至生態系統角度反映植被物候對氣候變化的響應。

遙感監測方法主要是將基于遙感植被指數時序數據(如NDVI、EVI等)識別出的物候期與氣候數據結合,分析氣候變化對物候的影響規律及影響程度[15,16]。遙感識別植被的物候期主要包含植被指數數據擬合和物候期求取。目前植被指數擬合方法主要有多項式擬合法、非對稱高斯法、分段邏輯斯蒂擬合法、Savitzky-Golay(S-G)濾波法等,物候期求取方法有閾值法、滑動平均法、曲率法、求導法等[17],應根據植被所在區域、植被類型等選取合適的方法來識別植被的物候期[18]以及進行其與氣候變化的相關性分析。例如,Piao等[19]利用1982-1999年NOAA/AVHRR NDVI數據,采用多項式擬合法對NDVI時序數據進行擬合,然后利用最大斜率閾值法識別中國溫帶植被的物候期,發現3-5月氣溫的升高導致返青期提前,8月中旬至10月氣溫的升高導致生長季結束日期延后;王植等[20]基于1982-2003年NOAA NDVI數據,采用Logistic函數擬合方法和曲率法研究了中國東部南北樣帶植被生長季起止日期,發現溫帶草原生長季起始和結束日期均呈提前趨勢,生長季起始日期與前一年冬季和當年春季溫度相關性較大,生長季結束日期與秋季降水相關;昝國盛等[21]利用1982-2006年NOAA/AVHRR NDVI數據,采用諧函數法和閾值法識別了呼倫貝爾草原植被的物候期,發現草原的生長季明顯延長,受春季溫度升高影響,返青期呈提前趨勢。遙感技術可以在群落或生態系統水平上監測植被物候對氣候變化的響應[22-24]。

目前,關于內蒙古羊草草原物候及其對氣候變化的響應研究大多數是基于地面觀測的物候數據進行的,少數基于遙感監測的研究,尚缺乏地面觀測的驗證,這使得遙感識別的物候期在數值和時空變化趨勢上的可靠性不強,從而影響植被物候變化趨勢的有效判斷和植被物候對氣候變化響應機理研究的深入開展。因此,本文基于地面觀測的物候數據、遙感監測NOAA NDVI數據和氣候數據,結合中國植被類型圖[2],選取覆蓋內蒙古羊草草原的實驗點,建立了經過農氣站點觀測結果驗證的物候期識別方法(NDVI曲線S-G濾波重建、滑動平均法求取物候期),獲得了1982-2013年內蒙古羊草草原的物候期,并結合各氣象站點及其對應年份的溫度、降水數據,分析了大范圍物候與主要氣象因子變化之間的作用關系。

1 數據來源

(1)地面物候數據。采用國家氣象局氣候中心提供的內蒙古地區各農氣站觀測的物候數據,選取了額爾古納市(50°15′N,120°11′E)、鄂溫克自治旗(49°09′N,119°45′E)、巴雅爾圖胡碩(45°04′N,120°20′E)和鑲黃旗(42°14′N,113°50′E)4個站點30年(1982-2011年)共113個樣本點(排除缺失值)的羊草返青期和黃枯期作為地面驗證數據。

(2)遙感數據。遙感數據來源于NOAA NDVI 3g數據(http://ecocast.arc.nasa.gov/data/pub/gimms/3g.v0/),其空間分辨率為8 km,時段為1982-2013年。該數據產品為15 d最大值合成,每年含有24個數值,32年的NDVI時序數據共含有768個數值。本實驗點是基于1∶100萬中國植被類型圖[2]選取的覆蓋內蒙古羊草草原且與農氣站和氣象站點鄰近的像元。

(3)氣象數據。氣象數據來源于中國氣象科學數據共享服務網(http://cdc.cma.gov.cn/),從其地面資料的中國地面氣候資料日值數據集中選取了覆蓋內蒙古羊草草原(或在農氣站附近)15個氣象站點的日均溫數據和降水數據,數據時段為1981-2013年(共33年)。

2 研究方法

2.1 技術路線

首先根據中國植被類型圖[2]中羊草草原的分布情況,選取內蒙古羊草草原的實驗點,從NOAA NDVI影像中提取各實驗點的1982-2013年NDVI時序數據;然后采用S-G濾波擬合法對每年NDVI時序數據進行擬合,將得到的NDVI時序曲線采用滑動平均法求取羊草草原的物候期(主要是返青期和黃枯期);最后利用地面物候觀測數據對求取的返青期和黃枯期進行驗證,以確定物候期識別方法的可行性。在此基礎上,結合日均溫數據和降水數據,進行遙感識別的物候期對氣候變化的響應分析。

2.2 NDVI曲線擬合去噪聲方法

Savitzky-Golay濾波法[25]是一種基于最小二乘原理的局域低階多項式擬合的數據平滑濾波方法,可在盡量保持原始信號規律的基礎上有效提高信噪比。其表達式為:

(1)

式中:yj、Yj分別表示第j個原始NDVI值和濾波后的NDVI值;Ci為局域內第i個樣本點的卷積系數;m為窗口寬度;共計m+1個原始樣本點參與第j個樣本點的低階(

2.3 物候期求取方法

滑動平均法[26,27]是利用原始植被指數曲線與其滑動平均曲線的交叉點判斷植被物候期。

Yt=(Xt+Xt-1+Xt-2+…+Xt-(w-1))/w

(2)

式中:Yt表示在t時刻的滑動平均值;Xt表示t時刻平滑后的NDVI值;w表示滑動平均的時間間隔(窗口寬度)。本文對比分析了180~290之間,間隔為10個DOY(Day Of Year)的12個窗口寬度,通過比較返青期、枯黃期識別結果的穩定性、遙感識別結果與地面觀測結果的一致性,確定返青期和枯黃期識別的最佳窗口寬度分別為280和190。

3 實驗結果及對比分析

3.1 物候期識別結果與驗證

圖1和圖2給出了4個氣象站30年(1982-2011年)遙感識別的返青期和黃枯期與地面物候觀測值的相關性分析結果。可以看出,基于S-G濾波擬合后采用滑動平均法識別返青期和黃枯期與地面觀測結果呈現出較顯著的正相關性,表明這種識別方法獲得的返青期和黃枯期變化趨勢與地面觀測結果一致,識別的返青期和黃枯期與地面觀測值之間的相關性均較高(相關系數(r)在0.43~0.63之間,顯著性水平(P)在0.005~0.05之間)。同時,識別的返青期和黃枯期與地面觀測結果在數值上差別主要集中在10~30 d,這可能是因為地面觀測的物候期是從植物個體尺度觀測的,而遙感監測的物候期是在群落甚至生態系統尺度監測的,監測方法的不同可能導致遙感監測結果較地面觀測結果有所滯后。

物候期識別方法的地面驗證結果表明,遙感物候期識別結果與農氣站觀測結果在物候變化趨勢上較一致。因此,本文采用S-G濾波擬合重建及滑動平均法識別羊草草原的返青期和黃枯期。

圖1 各氣象站遙感識別的返青期與地面觀測值的散點圖

圖2 各氣象站遙感識別的黃枯期與地面觀測值的散點圖

3.2 羊草草原物候及其生長季變化

在15個氣象站點附近選取1982-2013年NDVI數據,基于S-G濾波擬合方法重建NDVI曲線后,采用滑動平均法識別羊草草原的返青期和黃枯期。為分析相同年份各氣象站點物候期識別結果的統計分布情況,圖3-圖5給出了每年各站點物候期的四分位圖分析結果,為分析物候期隨年份的變化規律,對各站點遙感識別的物候期隨年份變化做了線性擬合分析,在圖3-圖5中增加了物候期隨年份變化的趨勢線及其顯著性評估結果。從圖3可以看出,1982-2013年各站點返青期主要集中在每年的第110-130天左右,同時發現近32年的返青期呈現出較顯著的提前趨勢,平均提前了2.2 d/10 a。從圖4可以看出,1982-2013年各站點黃枯期主要集中在每年的第255-275天左右,近32年的黃枯期變化趨勢不明顯。從圖5可以看出,1982-2013年各站點生長季長度主要為135~155 d左右,呈不顯著的延長趨勢,延長了1.4 d/10 a,這主要是由返青期提前導致的。

通過上述分析可知,同一年內,不同氣象站點的返青期、枯黃期、生長季長度均存在一定的差異,這與不同站點附近的遙感實驗點所處的地理環境及氣候條件有關,對于干旱、半干旱地區的草原植被而言,溫度、降水空間分布的差異可能是造成不同遙感實驗點物候期差異的重要原因。比較物候期的年際變化規律可知,1982-2013年返青期的提前趨勢較為顯著,黃枯期的變化趨勢相對不明顯,這表明返青期對氣候變化響應可能更為敏感,因此研究返青期與主要氣候因子的關系及其氣候響應機制,對于開展大面積物候變化預測、認識氣候變化對植被物候的影響等均具有重要參考價值。

圖3 各站點1982-2013年返青期變化趨勢

圖4 各站點1982-2013年黃枯期變化趨勢

圖5 各站點1982-2013年生長季長度變化趨勢

3.3 羊草草原物候對氣候變化的響應分析

表1給出了1982-2013年內蒙古15個氣象站點附近羊草草原的返青期、枯黃期、生長季長度變化量及其與溫度、降水的相關性分析結果。

1982-2013年大部分站點返青期提前7~10 d,少部分站點返青期提前1~4 d。返青期與春季溫度、冬春累積降水有明顯的相關性。返青期與春季氣溫,尤其是與當年4月、5月上旬平均氣溫呈負相關關系,與春季平均氣溫的相關系數多在-0.6~-0.2之間,與4月平均氣溫的相關系數多在-0.4~-0.1之間,與5月上旬平均氣溫的相關系數多在-0.7~-0.2之間。這種相關性的內在生長機制可能在于:羊草返青前需經歷一定程度的抗寒鍛煉及溫度驅動作用,內蒙古冬春溫度普遍較低,充足的抗寒鍛煉可以自然滿足,因而溫度驅動作用的強弱(即返青前春季氣溫的高低)便成為影響返青日期的重要氣象因素。降水對羊草草原的驅動效果也十分明顯,大部分站點返青期與前一年10月至當年4月的累積降水量呈負相關關系,相關系數多在-0.5~-0.2之間。這種相關作用的內在生長機制可能在于:在干旱/半干旱地區種子發育的休眠期[28],土壤墑情會對發育速度造成一定影響,在溫度滿足一定條件的前提下,充足的降水會對植被返青起到促進作用。

1982-2013年羊草草原的黃枯期也呈提前趨勢,但比返青期的提前幅度小,僅提前了0~4 d。黃枯期的變化主要與當年8月的平均氣溫呈負相關關系,相關系數在-0.5~-0.1之間,即黃枯前的氣溫越高,黃枯期越提前。

1982-2013年羊草草原的生長季長度整體呈延長趨勢,大部分站點延長了3~9 d。生長季長度的變化是返青和黃枯兩種發育過程綜合作用的結果,但主要是由返青期提前引起的,因為返青期比黃枯期提前的趨勢更加明顯。因此,生長季長度變化的主要影響因素與春季物候的影響因素相近,春季溫度及冬春累積降水作用較大,秋季溫度也有一定影響。

表1 1982-2013年內蒙古羊草草原物候對氣候變化的響應

4 結論與討論

本文基于NOAA NDVI數據和內蒙古部分氣象站點的羊草物候期地面觀測數據,在對原始植被指數進行去噪預處理[29]的基礎上,利用S-G濾波法進行了NDVI植被指數曲線的重建;采用滑動平均法進行了返青期和黃枯期的求取,并結合地面物候觀測結果進行了驗證,發現此法識別的物候期與地面觀測結果具有較好的一致性。然后,圍繞內蒙古羊草草原的15個氣象站點進行了遙感選點及物候期求取,并結合氣候數據對返青期、黃枯期和生長季長度的變化趨勢進行了分析,研究了物候期與氣候變化的相關性。結果表明:1982-2013年內蒙古羊草草原的返青期和黃枯期均呈提前趨勢,但返青期較黃枯期提前趨勢更加明顯,從而使生長季長度呈延長趨勢。返青期的提前與春季氣溫升高有關,尤其與當年4月和5月上旬氣溫的負相關作用明顯,與前一年10月至當年4月累積降水量增加也有顯著的負相關性;黃枯期的延后主要與當年8月氣溫的升高有關,與降水的相關性不大。

對于干旱/半干旱地區的草原植被而言,為更好地研究氣候對物候變化的影響,今后有必要考慮溫度和降水兩種主要氣象因素的聯合影響,同時考慮光照、土壤含水量等的綜合影響,最好是結合植物物候模型、地面物候觀測和遙感監測,從機理層面深入分析植物物候對氣候變化的響應。

[1] 張宏斌,唐華俊,楊桂霞,等.2000-2008年內蒙古草原MODIS NDVI時空特征變化[J].農業工程學報,2009,25(9):168-175.

[2] 中國科學院中國植被圖編輯委員會.中國植被圖集(1∶100萬)[M].北京:地質出版社,2007.

[3] 王旭陽,張顯峰,趙杰鵬,等.基于MODIS和TRMM數據的準格爾南緣植被凈初級生產力估算與分析[J].地理與地理信息科學,2011,27(2):21-25.

[4] DRAGONI D,SCHMID H P,WAYSON C A,et al.Evidence of increased net ecosystem productivity associated with a longer vegetated season in a deciduous forest in south-central Indiana,USA[J].Global Change Biology,2011,17(2):886-897.

[5] RICHARDSON A D,BLACK T A,CIAIS P,et al.Influence of spring and autumn phenological transitions on forest ecosystem productivity[J].Philosophical Transactions of the Royal Society B-Biological Sciences,2010,365(1555):3227-3246.

[6] WALTHER G R,POST E,CONVEY P,et al.Ecological responses to recent climate change[J].Nature,2002,416(6879):389-395.

[7] 楊柏娟,王思遠,常清,等.青藏高原植被凈初級生產力對物候變化的響應[J].地理與地理信息科學,2015,31(5):115-120.

[8] 顧潤源,周偉燦,白美蘭,等.氣候變化對內蒙古草原典型植物物候的影響[J].生態學報,2012,32(3):767-776.

[9] 楊曉華,越曉玲,娜日斯.內蒙古典型草原植物物候變化特征及其對氣候變化的響應[J].內蒙古草業,2010,22(3):51-56.

[10] 陳效逑,李倞.內蒙古草原羊草物候與氣象因子的關系[J].生態學報,2009,29(10):5280-5290.

[11] 張峰,周廣勝,王玉輝.內蒙古克氏針茅草原植物物候及其與氣候因子關系[J].植物生態學報,2008,32(6):1312-1322.

[12] 李榮平,周廣勝,王玉輝,等.羊草物候特征對氣候因子的響應[J].生態學雜志,2006,25(3):277-280.

[13] 翟佳,袁鳳輝,吳家兵.植物物候變化研究進展[J].生態學雜志,2015,34(11):3237-3243.

[14] MENZEL A,SPARKS T H,ESTRELLA N,et al.European phenological response to climate change matches the warming pattern[J].Global Change Biology,2006,12(10):1969-1976.

[15] 謝寶妮,秦占飛,王洋,等.基于遙感的黃土高原植被物候監測及其對氣候變化的響應[J].農業工程學報,2015(15):153-160.

[16] CONG N,WANG T,NAN H J,et al.Changes in satellite-derived spring vegetation green-up date and its linkage to climate in China from 1982 to 2010:A multimethod analysis[J].Global Change Biology,2013,19(3):881-891.

[17] HUDSON I L,KEATLEY M R.Phenological Research Methods for Environmental and Climate Change Analysis[M].Springer Dordrecht Heidelberg London New York:Springer,2010.

[18] 范德芹,趙學勝,朱文泉,等.植物物候遙感監測精度影響因素研究綜述[J].地理科學進展,2016,35(3):304-319.

[19] PIAO S L,FANG J Y,ZHOU L M,et al.Variations in satellite-derived phenology in China's temperate vegetation[J].Global Change Biology,2006,12(4):672-685.

[20] 王植.基于物候表征的中國東部南北樣帶上植被動態變化研究[D].北京:中國林業科學研究院,2008.

[21] 昝國盛,孫濤.呼倫貝爾草原植被覆蓋變化趨勢研究[J].林業資源管理,2011(1):44-48.

[22] XU H,TWINE T E,YANG X.Evaluating remotely sensed phenological metrics in a dynamic ecosystem model[J].Remote Sensing,2014,6(6):4660-4686.

[23] 陳效逑,王林海.遙感物候學研究進展[J].地理科學進展,2009,28(1):33-40.

[24] BADECK F W,BONDEAU A,B?TTCHER K,et al.Responses of spring phenology to climate change[J].New Phytologist,2004,162(2):295-309.

[25] CHEN J,J?NSSON P,TAMURA M,et al.A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter[J].Remote Sensing of Environment,2004,91(3-4):332-344.

[26] DUCHEMIN B,GOUBIER J,COURRIER G.Monitoring phenological key stages and cycle duration of temperate deciduous forest ecosystems with NOAA/AVHRR data[J].Remote Sensing of Environment,1999,67(1):68-82.

[27] REED B C,BROWN J F,VANDERZEE D,et al.Measuring phenological variability from satellite imagery[J].Journal of Vegetation Science,1994,5(5):703-714.

[28] CESARACCIO C,SPANO D,SNYDER R L,et al.Chilling and forcing model to predict bud-burst of crop and forest species[J].Agricultural and Forest Meteorology,2004,126(1-2):1-13.

[29] 范德芹,朱文泉,潘耀忠,等.基于狄克松檢驗的NDVI時序數據噪聲檢測及其在數據重建中的應用[J].遙感學報,2013,17(5):1158-1174.

關于防范虛假投稿網站的公告

《地理與地理信息科學》是地理類綜合性學術期刊(雙月刊),全國中文核心期刊、中國科技核心期刊、中國科學引文數據庫(CSCD)核心期刊,國內外公開發行。它的任務是:堅持面向經濟建設,面向科研開發,重點反映地理學與地理信息科學的新理論、新技術、新方法、新成果。于2016年正式啟用作者在線投稿系統,投稿網址為:http://dlgt.cbpt.cnki.net,歡迎廣大作者不吝賜稿。

近期,有人以我社的名義,通過虛假網站或郵件等方式,接收投稿并收取費用,此舉對《地理與地理信息科學》雜志社及作者造成了不良影響。

我社特此聲明如下:我社不會以任何形式要求作者向個人賬戶匯款,請注意甄別。如遇到此類情況,歡迎廣大作者向我社提供相關線索!《地理與地理信息科學》編輯部郵箱為:dlxxkx@vip.163.com;聯系電話:0311-86054904;通訊地址:石家莊市槐安東路116號/育才街252號北樓101室。

《地理與地理信息科學》雜志社

2016年11月15日

Phenology ofLeymusChinensisSteppe in Inner Mongolia and Its Response to Climate Changes

FAN De-qin1,ZHAO Xue-sheng1,ZHENG Zhou-tao2

(1.CollegeofGeosciencesandSurveyEngineering,ChinaUniversityofMining&Technology(Beijing),Beijing100083;2.CollegeofResourcesScienceandTechnology,BeijingNormalUniversity,Beijing100875,China)

Along with the global climate change,plant phenology has also changed.In order to analyze the response of phenology forLeymuschinensissteppe in Inner Mongolia to climate effectively,we first used the Savitzky-Golay function to reconstruct the National Oceanic and Atmospheric Administration Normalized Difference Vegetation Index (NDVI) time series for the selected test pixels,then used moving average methods to derive the phenology forLeymuschinensissteppe.The performance of the methods for green-up date and the end of season date identification were verified by the ground observation samples.The phenology dates identified by the moving average method using the reconstructed NDVI time series,were well consistent with the observed ground phenology data.The correlations coefficient (r) of the phenology derived from remote sensing and the ground observation were between 0.43 and 0.63,while the significant level (P) were between 0.005 and 0.05.The differences of the phenology derived from remote sensing and the ground observation were between 10 and 30 days.On this basis,we acquired phenology stages (e.g.green-up dates and end of season dates) and the length of the growing season at each weather station,and analyzed their response to climate change combing with climate data (temperature and precipitation data) from 1982 to 2013.The result showed that the mean green-up dates forLeymuschinensissteppe in Inner Mongolia were mainly occurred between DOY (Day Of Year) 110 and 130,and the green-up dates were advanced at a rate of 2.2 days per decade during 1982-2013.The mean end of season dates forLeymuschinensissteppe in Inner Mongolia were mainly occurred between DOY 255 and 275,which were slightly advanced during 1982-2013.The trend of green-up dates forLeymuschinensissteppe in Inner Mongolia advanced,which was related with the rising spring temperature and the increasing accumulated precipitation between October of the last year and April of this year.The trend of end of season dates also advanced,which was related with the rising temperature of this year,but not related with the precipitation.Therefore,the length of growing season was delayed mainly due to the advanced trend of green-up dates was more obviously than the end of season dates.

Leymuschinensis;steppe;NDVI;phenology;remote sensing;climate

2016-06-27;

2016-08-26

國家自然科學基金項目(41171306、41601456)

范德芹(1982-),女,博士后,主要從事植被與生態遙感研究。*通訊作者E-mail: zxs@cumtb.edu.cn

10.3969/j.issn.1672-0504.2016.06.014

TP79

A

1672-0504(2016)06-0081-06

主站蜘蛛池模板: 国产一级片网址| 国产一级二级三级毛片| 波多野结衣在线se| 久久久久九九精品影院 | 91久久大香线蕉| 天天爽免费视频| 国产午夜福利在线小视频| 男人天堂亚洲天堂| 国产亚洲欧美在线专区| 国产精品自在在线午夜| 亚洲熟妇AV日韩熟妇在线| 美女无遮挡被啪啪到高潮免费| 亚洲国产理论片在线播放| 暴力调教一区二区三区| 亚洲日韩高清在线亚洲专区| 丰满人妻中出白浆| 88国产经典欧美一区二区三区| 五月天福利视频| 九九热这里只有国产精品| 亚洲制服丝袜第一页| 色视频国产| 亚洲国产成人麻豆精品| 欧美国产日韩在线播放| 一本大道视频精品人妻 | 91精品伊人久久大香线蕉| 国内精品免费| 一级成人a做片免费| 国产三级a| 欧美综合激情| 1769国产精品免费视频| 亚洲欧美国产高清va在线播放| 免费一看一级毛片| 国产精品亚欧美一区二区| 国产一区二区三区精品欧美日韩| 综合人妻久久一区二区精品 | www.亚洲国产| 九九九久久国产精品| 久久精品中文字幕少妇| 毛片网站观看| 91精品福利自产拍在线观看| 成人在线天堂| 国产嫖妓91东北老熟女久久一| 国产成人综合日韩精品无码首页 | 2020国产免费久久精品99| 久久综合国产乱子免费| 成人国产一区二区三区| 国产熟女一级毛片| 特级欧美视频aaaaaa| 2020精品极品国产色在线观看 | 国产一区二区三区精品久久呦| 成人综合在线观看| 狠狠躁天天躁夜夜躁婷婷| 欧美一级视频免费| 久久婷婷五月综合色一区二区| 女人天堂av免费| 久久综合丝袜日本网| 国产无人区一区二区三区| 久久精品电影| 国产成人精品一区二区免费看京| 高清欧美性猛交XXXX黑人猛交 | 69视频国产| 伊人国产无码高清视频| 久热中文字幕在线| 久久毛片基地| 极品国产在线| 免费人成网站在线观看欧美| 亚洲日韩精品欧美中文字幕| 中文字幕2区| 一级毛片免费观看久| 欧美不卡二区| 久久香蕉国产线看观看精品蕉| 精品国产欧美精品v| 老汉色老汉首页a亚洲| 欧美三级不卡在线观看视频| 五月婷婷综合色| 久久精品中文字幕免费| 国产电话自拍伊人| 狼友视频一区二区三区| 四虎成人精品| 黄色一级视频欧美| 国产亚洲精品97AA片在线播放| 欧美va亚洲va香蕉在线|