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

基于AMSR-E數據的微波植被指數與MODIS植被指數關系研究

2012-09-18 06:32:55陳思宇于惠馮琦勝呂志邦梁天剛
草業科學 2012年3期
關鍵詞:研究

陳思宇,于惠,馮琦勝,呂志邦,梁天剛

(草地農業生態系統國家重點實驗室蘭州大學草地農業科技學院,甘肅蘭州730020)

基于AMSR-E數據的微波植被指數與MODIS植被指數關系研究

陳思宇,于惠,馮琦勝,呂志邦,梁天剛

(草地農業生態系統國家重點實驗室蘭州大學草地農業科技學院,甘肅蘭州730020)

利用青南牧區2007-2010年的AMSR-E亮溫數據計算了相應的微波植被指數(Microwave Vegetation Index,MVI),對MVI的月季動態變化特征進行了分析,同時結合相同時間序列的MODIS NDVI和EVI數據,對比分析了MVI和MODIS植被指數之間的相關關系,篩選出NDVI反演模型,并對模型的精度進行了評價。結果表明,MVI值隨著植被的生長而降低;MVI與NDVI、EVI均有顯著的線性負相關。其中,升軌低頻MVI與NDVI的相關性最好,相關系數為0.58(P<0.001);MVI與MODIS植被指數之間的最優模型為NDVI=-0.85×MVI+0.84;利用最優模型將反演的NDVI與MODIS NDVI進行比較,兩者差異較小,說明這一模型能較好地反映2種植被指數的關系。

微波植被指數;MODIS植被指數;相關性;反演模型

植被指數作為評價植被覆蓋度、生長活力和生物量等植被信息的重要手段,已廣泛應用于許多研究領域[1-3]。據統計,經過40多年的研究,在光學遙感領域已發展了40多種植被指數。但是光學遙感容易受到植被本身、環境條件和大氣狀況等多種因素的影響,在時間和空間上均具有一定的局限性。微波遙感不受太陽照射、大氣、云層、降水等因素的限制,具有全天時和全天候工作的能力,對植被木質部生物量、植被類型以及土壤含水量等具有非常強的敏感性,可以探測到相對較厚的植被層,這是其他遙感手段所不能得到的[4-9]。因此,利用微波傳感器不同頻率的亮溫數據,獲取陸表植被參數信息,在彌補光學植被指數的缺陷方面具有重要的科學意義[10]。

盡管微波遙感相對于光學遙感存在諸多優勢,但多年來在植被探測方面的研究應用較少[11]。Paloscia和Pampalonil[12]通過10GHz和36GHz水平以及垂直4個通道的亮溫數據對農作物生長進行了監測,結果表明,微波植被指數(Microwave Vegetation Index,MVI)對植被類型以及植被覆蓋下的土壤水分非常敏感,頻率10GHz的微波極化指數(Microwave Polarization Index,MPI)與葉面積指數(Leaf Area Index,LAI)之間存在指數關系。比較分析MVI與NDVI在灌木、農田和草原3種土地覆蓋類型下的相關關系,發現MVI能夠更好地反映植物木質部生物量、結構、含水量以及鮮生物量等信息[13]。呂京國等[14]在全國范圍內選取5種典型地表類型分析微波植被指數與NDVI的關系,結果表明,兩者呈負相關關系,相關系數為-0.46。此外,Choudhury和Tucker[15-16]、Ulaby等[17]、Jackson和Schmugge[18]等,也在主動微波指數與植被監測方面做了一些相關研究。可以看出,目前對于微波植被指數的研究,主要集中在微波指數與光學指數的定性分析方面,而對于它們之間的關系以及關系的數學表達等方面的工作還較少。因此,為了更好地利用微波植被指數監測植被動態變化,需要對其進行更加深入的研究。

本研究在前人工作的基礎上,以青南牧區為例,利用同一時間序列的MVI和MODIS植被指數數據,研究MVI與NDVI和EVI之間的相關性,構建NDVI反演模型,并對其精度進行綜合評價,以期為更加科學有效地監測植被的動態變化提供科學依據。

1 數據與方法

1.1 研究區概況 青南牧區地處青藏高原腹地,位于31°38′~36°20′N,89°31′~102°14′E,是我國黃河、長江、瀾滄江三大水系發源地。青南牧區土地總面積3.56×107hm2,約占青海省土地總面積的50.4%。平均海拔4 000m以上,年平均氣溫為-5.6~4.9℃,≥0℃的年積溫一般不超過1 500℃·d,不少地區低于1 000℃·d,年降水量基本在391.7~764.0mm。草地是青南牧區土地資源的主體,現有天然草原面積2.10×107hm2,占青南牧區土地總面積的59%,其中可利用草地面積1.79×107hm2。草地類型主要有高寒草甸類和高寒干草原類,分別約占草地總面積的76%和23%[19-20]。

1.2 研究數據 主要包括2類遙感資料:1)高級微波掃描輻射計(Advanced Microwave Scanning Radiometer-Earth Observing System,AMSR-E)數據:AMSR-E是搭載在美國EOS-AQUA衛星上的圓錐掃描方式微波成像儀,由NASA(National Aeronautics and Space Administration,NASA)于2002年5月4日發射升空。AMSR-E傳感器的工作頻段有6個(6.9~89.0GHz),有水平和垂直2種極化模式,共12個通道。此外,每天根據過境時間的不同又有升軌和降軌2次數字圖像。本研究使用6GHz、10GHz和18GHz水平與垂直極化的6個通道的亮溫數據,覆蓋范圍為北半球,投影格式為北半球可擴展的等面積地球格網(EASE-Grid_north),像元大小25km,時間序列為2007-2009年以及2010年4-10月的每日升軌和降軌數據,共有亮度溫度圖像16 440幅[21]。2)MODIS月最大合成植被指數產品(MOD13A3):由對地觀測系統數據共享平臺(EOS Data Gateway,EDG)提供,空間分辨率1km。青南牧區每月包括兩幅數字圖像,編號分別為H25V05和H26V05,時間序列2007-2009年1-12月及2010年4-10月,共計36景EVI圖像和43景NDVI圖像。

1.3 AMSR-E和MODIS數據處理方法AMSR-E亮度溫度數據(Tb,Brightness Temperature)和MOD13A3數據的處理流程如圖1所示,具體步驟如下:1)AMSR-E Tb數據處理。首先,給解壓的原文件加.BSQ擴展名,并且對每個BSQ文件建立相應的頭文件;其次,在ArcGIS軟件的ArcMap模塊中將文件轉換為GRID格式,并將投影定義為EASE-Grid_north,利用ArcInfo工作站軟件將投影轉換為Albers Krasovsky,將格網大小重采樣到25 000m;最后分別提取逐日升軌及降軌6個通道的亮溫值。2)MOD13A3月植被指數產品處理。利用MODIS數據重投影工具(MODIS Reprojection Tools,MRT)將HDF文件轉為TIF格式,在ArcMap中轉換為GRID格式,使用Arc-Info工作站將投影轉為Albers Krasovsky,重采樣到1 000m。

1.4 微波植被指數算法 裸露地表情況下,不同頻率輻射率之間呈很強的線性關系。根據這一特點,利用地表輻射模型,可推導出微波植被指數[11]。

圖1 AMSR-E Tb數據和MOD13A3數據處理流程圖Fig.1 Data processing for AMSR-E Tb and MOD13A3

式中,TBv(f2)、TBh(f2)分別表示頻率f2的垂直和水平極化的亮溫值,TBv(f1)、TBh(f1)分別表示頻率f1的垂直和水平極化亮溫值。f1、f2是亮溫數據中的相鄰頻率(f2>f1)。

根據式(1),選取6、10和18GHz水平和垂直極化6個通道的亮溫數據,計算6和10GHz及10和18GHz 2組相鄰頻率梯度下的MVI值。此外,因每日有升軌和降軌2次數字圖像。因此,本研究共計算升軌低頻、升軌高頻、降軌低頻、降軌高頻4種微波植被指數,分別用MVI(6,10)A、MVI(10,18)A、MVI(6,10)D、MVI(10,18)D來表示。

1.5 數據篩選方法 在無大氣因素的影響下,除積雪覆蓋的區域外高頻微波信號通常強于低頻微波信號,但射頻干擾(Radio-Frequency Interference,RFI)導致了不規則的頻率間梯度變化,這種不規則的梯度差通常導致MVI值異常。為了防止RFI信號和積雪對研究區域數據的影響,需要對部分受到射頻干擾的微波數據進行篩選。本研究數據篩選的方法為:1)水平極化通道的微波亮溫數據在受RFI干擾后,其值將大于垂直極化通道的微波亮溫數據,據此剔除水平極化亮溫值大于垂直極化亮溫值的異常值;2)受到RFI的干擾后,高頻和低頻亮溫數據間將呈負向梯度,為了避免這種干擾剔除高頻和低頻梯度差小于等于-5的異常值;3)由于射頻干擾的影響,MVI值將會超出其正常范圍,因此刪除MVI值小于0及大于1的數據[22]。

研究區升軌和降軌MVI數據的原始樣本數各有1 353 399個,按上述所列條件篩選后,升軌樣本數為660 826個,降軌樣本數為383 241個。由此可見,RFI對研究區被動微波數據有嚴重的影響。

1.6 月MVI數據的合成 對微波植被指數而言,當植被密度增加時,微波穿透植被層的能力減弱,MVI值減小。不同的植被類型,由于其枝干形狀、葉片朝向等的不同,其MVI值也是不同的。同時,同一植被類型在物候循環中的不同階段,其生命力特征及表現形式也不同。而NDVI與EVI值會隨著地表植被密度的增加而增大[13]。

為了對植被指數進行比較分析,必須將植被指數數據在時間和空間尺度上進行匹配。本項研究采用最小值合成法將計算得到的每日MVI數據合成月MVI數據,具體方法如下:根據式(1)計算出每天的MVI,然后用最小值法將每天的MVI合成為月最小MVI數據。利用ArcGIS對MODIS月植被指數NDVI和EVI進行區域統計,分別計算出在25 km格網中NDVI和EVI的平均值。在此基礎上,分析比較MODIS的2種植被指數同MVI之間的關系。

2 結果與討論

2.1 MVI的月季變化特征 按照MVI的不同,曲線分為升軌低頻、升軌高頻、降軌低頻和降軌高頻4種。青南牧區2007-2009年月MVI平均值的動態曲線總體趨勢相同,且變化平緩(圖2)。

圖2 微波植被指數(MVI)月動態變化圖Fig.2 Dynamic of monthly MVI

隨著牧草在4月開始返青,曲線逐漸呈下降趨勢,在盛草時期(7-9月)植被生長達全年的最大值,相應的MVI曲線出現低谷,10月以后草地開始枯黃,牧草進入休眠期,MVI曲線開始上升(圖2)。高頻MVI值在全年內均小于低頻MVI值,這是因為高頻微波穿過植被層時,散射現象較低頻微波的嚴重,從而到達傳感器的能量較低頻微波的弱。

2.2 微波植被指數和MODIS植被指數的相關分析 從研究區3年MVI與NDVI和EVI之間的散點圖可以看出(圖3),MVI與MODIS植被指數之間有明顯的線性關系。統計分析的結果表明,MVI與MODIS植被指數均呈負相關關系(P<0.001)。當NDVI和EVI較大時,MVI較小;相反,當NDVI和EVI很小(即植被狀況很差或為裸地)時,MVI很大。無論升軌還是降軌數據計算的MVI,其低頻微波植被指數分布相對集中,而高頻指數較為分散。這說明低頻微波植被指數與MODIS植被指數的相關關系優于高頻微波植被指數。在升軌時期植被受水分因素的影響較小,這可能是升軌MVI優于降軌的主要原因。

圖3 2007-2009年NDVI、EVI與MVI相關性散點圖Fig.3 Correlation between NDVI,EVI,and MVI from 2007to 2009

MVI與NDVI和EVI的相關系數差別不大,在0.505~0.588波動(圖3)。低頻的MVI與MODIS植被指數的相關性優于高頻的MVI,尤其是升軌的MVI,其低頻優于高頻的趨勢更加明顯。升軌低頻的MVI與NDVI的相關關系最好,相關系數達0.588,而升軌高頻的MVI與NDVI的相關關系最差,相關性系數為0.505,其余的相關系數介于二者之間。由于NDVI季節性變化較明顯,主要反映植被表層的覆蓋狀況,而MVI可反映植被結構等方面的特征信息。因此,MVI隨著植被生長周期性變化小,而NDVI隨植被生長的周期性變化較大,從而導致兩者的相關關系受到影響,造成相關系數偏低[13]。另外,AMSR-E數據空間分辨率為25km,遠低于MODIS 1km的分辨率,這也是造成兩者相關性較低的另一重要原因。

2.3 精度分析 利用青南牧區2010年4-10月的亮溫數據計算低頻升軌的MVI值,用低頻升軌的MVI值和相應的統計模型(圖3)反演出研究區NDVI值,并將反演的NDVI值與MODIS的NDVI值進行比較,驗證模型精度。在4-10月中,6月的均方根誤差(Root Mean Square Error,RMSE)最小,是實測值與預測值之間的離散程度最低的月份(表1)。實測NDVI的月最大值的波動范圍較反演的NDVI值大。匯總7個月的實測值和預測值發現,4、5月實測的NDVI最大值和平均值均小于預測NDVI的最大值和平均值,7-10月的實測值最大值和平均值均大于預測值的最大值和平均值。總體來說,利用反演模型計算的NDVI值低于實際的NDVI值,均方根誤差為0.158。青南牧區地形復雜,海拔梯度較大。9月,受地形的影響研究區內牧草長勢不均一,在高海拔地區牧草已經進入枯黃期,而低海拔地區的牧草可能還處于生長旺季。NDVI對植物體內葉綠體變化敏感,因而NDVI值在整個研究區內的變化范圍較大。MVI反演出的NDVI值相對穩定,變化趨勢平緩。這可能是導致9月二者之間的離散程度達到最大的一個重要原因,均方根誤差達0.235。從NDVI實測值和模擬值之間的關系散點圖可以看出,NDVI的實測值與模擬值呈現較好的正相關關系(圖4),其相關系數為0.653。

表1 MODIS NDVI實測值與預測值Table 1 Observed and pridicted value of MODIS NDVI

圖4 MODIS NDVI數據值與預測值之間的關系Fig.4 Relationship between observed and predicted value of MODIS NDVI

利用研究區植被生長狀況最好時期8月的MODIS NDVI影像及NDVI模擬圖像,合成實測NDVI和模擬NDVI的差值圖(圖5)。研究區西部屬東昆侖山地區,該區平均海拔較高,草地類型主要為山地草甸類和高寒草原類。這2類草地的共同特點是植被覆蓋度較低,在此區域實測NDVI與預測NDVI的差值主要介于-0.20~0.00;在中部地區,海拔差異大,地形起伏劇烈,地表植被類型復雜且草地比較分散,二者間的差值參差不齊,主要集中在-0.20~0.00、0.00~0.20以及0.20~0.40,個別地區的差值在0.40~0.60;在東部地區,二者差值主要分布于0.00~0.20和0.20~0.40。東部地區海拔較低且草地類型均一,主要為高寒草甸類,此類草地的覆蓋度為80%~90%,植被狀況較好,其預測的NDVI值偏低。從NDVI實測值與預測值之間的差值頻數分布狀況(圖6)可以看出,整個研究區實測NDVI和預測NDVI的差值主要分布在0.00~0.20以及0.20~0.40。總體來說,模型預測的NDVI值比較接近MODIS NDVI數據值。

圖5 青南牧區實測與模擬NDVI差值圖Fig.5 Differenc in observed and simulated value of NDVI in pastoral areas of southern Qinghai

圖6 實測與模擬NDVI差值分布直方圖Fig.6 Histogram of the difference in observed and estimated values of NDVI

3 結論與展望

本研究表明,微波植被指數與MODIS植被指數之間有較好的相關性,隨著微波植被指數的減小,NDVI和EVI逐漸增大。低頻MVI與NDVI、EVI的相關關系優于高頻的MVI,高頻微波在穿透植被層時散射現象比低頻微波嚴重。利用2010年4-10月的MVI和MODIS NDVI數據,對最優反演模型進行了精度驗證,結果表明模型能較好地反映MODIS植被指數和微波植被指數之間的關系。

被動微波數據時間分辨率高且不受天氣狀況的影響,是植被監測的重要手段。但其空間分辨率較低(25km),混合像元的問題比較嚴重,這一問題是當前研究的難點所在。因此,今后的研究應致力于如何有效結合多方面的資料,消除混合像元對微波數據的影響,從而更加準確地監測植被生長狀況,獲取更多的植被信息。

[1]孫紅雨,王長耀,牛錚,等.中國地表植被覆蓋變化及其與氣候因子關系——基于NOAA時間序列數據集[J].遙感學報,1998,2(3):205-210.

[2]王建偉,陳功.草地植被指數及生物量的估測[J].云南農業大學學報,2006,21(3):372-375.

[3]李海亮,趙軍.草地遙感估產的原理和方法[J].草業科學,2009,26(3):34-38.

[4]張凱,郭鈮,王潤元,等.甘南草地地上生物量的高光譜遙感估算研究[J].草業科學,2009,26(11):44-50.

[5]王鶯,夏文韜,梁天剛,等.基于MODIS植被指數的甘南草地凈初級生產力時空變化研究[J].草業學報,2010,19(1):201-210.

[6]田慶久,閔祥軍.遙感植被指數研究進展[J].地球科學進展,1998,13(4):227-233.

[7]陳亮,杜金陽.L波段多角度微波植被指數研究[J].理論研究,2010(3):13-16.

[8]米兆榮,張耀生,趙新泉,等.NDVI和EVI在高寒草地牧草鮮質量估算和植被動態監測中的比較[J].草業科學,2010,27(6):13-19.

[9]毛克彪,唐華俊,周清波,等.AMSR-E微波極化指數與MODIS植被指數關系研究[J].國土資源遙感,2007(1):27-31.

[10]張鐘軍,孫國清,張立新,等.被動微波遙感中一種基于輻射傳輸理論的植被層模型[J].北京師范大學學報(自然科學版),2003,39(5):694-700.

[11]高峰,車濤,王介民,等.被動微波遙感指數極其應用[J].遙感技術與應用,2005,20(6):551-557.

[12]Paloscia S,Pampalonil P.Microwave vegetation indexes for detecting biomass and water conditions of agricultural crops[J].Remote Sense of Environment,1992,40:15-26.

[13]Shi J C,Jackson T,Tao J,etal.Microwave vegetation indices for short vegetation covers from satellite passive microwave sensor AMSR-E[J].Remote Sense of Environment,2008,112:4285-4300.

[14]呂京國,張小詠,蔣玲梅,等.微波穿透指數MVI與光學植被指數NDVI的關系探討[J].遙感應用,2009(6):39-42.

[15]Choudhury B J.Relationship between vegetation indices,radiation absorption,and net photosynthesis evaluated by a sensitivity analysis[J].Remote Sense of Environment,1987,22:209-223.

[16]Choudhury B J,Tucker C J.Monitoring global vegetation using Nimbus-7 37GHz data:Some empirical relations[J].International Journal of Remote Sensing,1987,8(7):1085-1090.

[17]Ulaby F T,Razani M,Dobson M C.Effects of vegetation cover on the micro-wave radiometric sensitivity to soil moisture[J].IEEE Transations on Geoscience and Remote Sensing,1983,21:51-61.

[18]Jackson T J,Schmugge T J.Vegetation effects on the microwave emission from soils[J].Remote Sense of Environment,1991,36:203-212.

[19]郝云晴,于明勝.青南牧區草地生態牧業建設的思考及對策[J].青海畜牧獸醫雜志,2002,32(4):33-35.

[20]趙小娟.青南牧區生態畜牧業發展現狀及對策[J].青海畜牧獸醫雜志,2008,38(5):34-36.

[21]于惠,馮琦勝,張學通,等.基于AMSR-E信息的北疆牧區雪深遙感監測模型方法初探[J].草業學報,2009,18(4):210-216.

[22]Li L,Njoku E,Im E,etal.A Preliminary Survey of Radio-frequency Interference over the U.S.in Aqua AMSR-E Data[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(2):380-390.

Relationship between microwave vegetation index based on AMSR-E data and vegetation index based on MODIS data

CHEN Si-yu,YU Hui,FENG Qi-sheng,LV Zhi-bang,LIANG Tian-gang
(State Key Laboratory of Grassland Agro-ecosystems,College of Pastoral Agriculture Science and Technology,Lanzhou University,Lanzhou 730020,China)

Microwave Vegetation Indices(MVIs)in the pastoral area of the south of Qinghai Province were calculated by using AMSR-E brightness temperature data from 2007to 2010and its features of monthly MVI changes were discussed,and then the relationship between microwave vegetation index based on AMSR-E data and vegetation index based on MODIS data for same time series was determined in this study.This study showed that the MVI value reduced as the vegetation plant grew and the MVIs was strongly negative correlations with NDVI and EVI,in which the correlation between MVIs and NDVI was the best with the correlation coefficient of 0.58.The optimal model between MVIs and MODIS vegetation indices wasNDVI=-0.85×MVI+0.84,and the comparison result showed that the difference between MODIS NDVI and simulated NDVI from the optimal model was little.This study suggested that the optimal model reflected the relationship of the two vegetation indices.

Microwave Vegetation Index;MODIS vegetation index;correlation;simulating model

LIANG Tian-gang E-mail:tgliang@lzu.edu.cn

TP79;Q94

A

1001-0629(2012)03-0377-07

2011-03-31 接受日期:2011-07-11

教育部高等學校科技創新工程重大項目培育資金項目(708089);國家科技支撐計劃項目(2009BAC53B01);國家高技術研究發展專項(2007AA10Z232)

陳思宇(1987-),女,甘肅定西人,在讀碩士生,研究方向為草地遙感與地理信息系統。E-mail:chensy_10@lzu.edu.cn

梁天剛 E-mail:tgliang@lzu.edu.cn

猜你喜歡
研究
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
關于遼朝“一國兩制”研究的回顧與思考
EMA伺服控制系統研究
基于聲、光、磁、觸摸多功能控制的研究
電子制作(2018年11期)2018-08-04 03:26:04
新版C-NCAP側面碰撞假人損傷研究
關于反傾銷會計研究的思考
焊接膜層脫落的攻關研究
電子制作(2017年23期)2017-02-02 07:17:19
主站蜘蛛池模板: 国产日本欧美亚洲精品视| 欧美另类第一页| 亚洲天堂久久| 久久综合干| 亚洲精品另类| 免费国产好深啊好涨好硬视频| 三区在线视频| 欧美日韩v| 奇米精品一区二区三区在线观看| 91精品免费高清在线| 91人妻在线视频| 中文国产成人精品久久| 中国精品久久| 青青青国产视频| 久久夜夜视频| 热九九精品| 成人在线观看不卡| 精品无码人妻一区二区| 成人午夜视频网站| 亚洲综合专区| 成人看片欧美一区二区| 色播五月婷婷| 夜夜操国产| 欧美翘臀一区二区三区| 青青青伊人色综合久久| 国产超碰在线观看| 国产va免费精品观看| 精品国产免费第一区二区三区日韩| 亚洲国产成人自拍| 99re精彩视频| 97在线碰| 亚洲精品视频网| 亚洲五月激情网| 日本人又色又爽的视频| 免费三A级毛片视频| 综合五月天网| 亚洲自拍另类| 好吊妞欧美视频免费| 国产91特黄特色A级毛片| 国产一级毛片高清完整视频版| 精品福利网| 国产精品视频a| 国产成人狂喷潮在线观看2345| av在线人妻熟妇| 99久久精品免费视频| 国产一二三区在线| 92午夜福利影院一区二区三区| 亚洲国产精品一区二区第一页免| 日本三级黄在线观看| 亚洲人成网址| 狠狠色丁香婷婷| 天堂网亚洲系列亚洲系列| 爱爱影院18禁免费| 一区二区理伦视频| 影音先锋亚洲无码| 免费欧美一级| 伊人成人在线视频| 91九色国产porny| 无码内射中文字幕岛国片| 一本色道久久88综合日韩精品| 免费观看精品视频999| 国产美女精品人人做人人爽| 亚洲第一页在线观看| 国产在线观看第二页| 欧美日韩国产精品综合| 中文字幕乱码中文乱码51精品| 亚欧成人无码AV在线播放| 中字无码av在线电影| 91视频免费观看网站| 日日拍夜夜嗷嗷叫国产| 欧美激情伊人| 久久久久国产精品嫩草影院| 欧美色视频日本| 亚洲精品无码高潮喷水A| 91国内在线观看| 美女视频黄又黄又免费高清| 午夜少妇精品视频小电影| 色悠久久久久久久综合网伊人| 久热这里只有精品6| 国模私拍一区二区| 国产系列在线| 精品黑人一区二区三区|