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

利用溫濕度指數提高紅樹林遙感識別精度

2012-12-27 07:30:52張雪紅田慶久
自然資源遙感 2012年3期
關鍵詞:紅樹林分類特征

張雪紅,田慶久

(1.南京信息工程大學遙感學院,南京 210044;2.南京大學國際地球系統科學研究所,南京 210093)

利用溫濕度指數提高紅樹林遙感識別精度

張雪紅1,田慶久2

(1.南京信息工程大學遙感學院,南京 210044;2.南京大學國際地球系統科學研究所,南京 210093)

針對使用TM圖像反射波段信息難以將紅樹林與陸地植被、尤其是與水體—植被混合像元有效區分的問題,結合不同潮位的TM圖像,基于反射波段信息,引入TM6熱紅外波段信息,提出了溫濕度指數(temperature-moisture index,TMI)。分析結果表明,綜合潮位信息、熱紅外波段信息及溫濕度指數能顯著提高紅樹林與其他地物之間的可分性。采用光譜角度制圖(spectral angle mapping,SAM)監督分類法對紅樹林進行分類識別,較之其他研究者所采用的分類特征,熱紅外波段信息及溫濕度指數能使紅樹林分類精度明顯提高(Kappa系數提高了0.14,錯分率降低了19.9%),說明利用潮位信息、熱紅外波段信息及溫濕度指數可以提高紅樹林的遙感識別精度。

紅樹林;潮位;TM;熱紅外;溫濕度指數

0 引言

紅樹林是生長在熱帶、亞熱帶沿海潮間帶灘涂上特有的木本植物群落,屬常綠闊葉林,主要分布于淤泥深厚的海灣或河口鹽漬土壤上。紅樹林具有促淤固灘、防浪護堤、保持生物多樣性和凈化環境污染等作用[1-2]。紅樹林的特殊生長環境不利于大面積的野外實地調查,而需借助遙感技術進行快速、大面積地監測和調查。關于紅樹林的遙感識別與監測,已有研究的關注重點主要包括2個方面:①紅樹林分布的面積與變化,側重于從宏觀尺度上監測紅樹林的分布范圍及動態變化,故此類研究多采用Landsat MSS[3],Landsat TM(ETM+)[3-7],SPOT[5,8-10]和ASTER[11]等中分辨率衛星遙感圖像數據;②紅樹林的種類及其群落分布規律,以及生物量、葉面積指數等生物物理參數,主要在微觀尺度上研究紅樹林生態群落分布及演變規律,故此類研究多采用航空遙感圖像(如 CASI圖像[6],HyMap 圖像[12]),IKONOS[13-15]以及 QuickBird[13,16]等高分辨率或高光譜圖像數據。

在紅樹林的識別研究中,主要采用 NDVI[17],TM3/TM5[6],(TM5 - TM7)/(TM5+TM7)[4]和TM5/TM4[6]等植被指數或波段比值組合來區分紅樹林與非紅樹林。部分研究中結合了DEM和海岸線等輔助地學數據,以提高紅樹林的識別精度[2]。由于紅樹林屬于濕地環境中的植物,使得紅樹林與水體—植被混合像元的反射光譜極其相似,因此紅樹林與陸地植被、尤其是與水體—植被混合像元之間存在嚴重的“同譜異物”現象。僅僅依靠遙感圖像的反射波段信息常常難以完全消除這種“同譜異物”效應,嚴重制約著紅樹林的遙感識別精度。

針對紅樹林與水體—植被混合像元難以區分的問題,本文從它們的“同譜異物”效應出發,基于紅樹林的濱海環境特點,在利用遙感圖像反射波段信息的基礎上,引入能反映它們之間差別的熱紅外波段信息,并提出溫濕度指數(temperature-moisture index,TMI),從而提高了紅樹林與水體—植被混合像元的可分性,使“同譜異物”轉變成“異譜異物”,進而提高紅樹林的遙感識別精度。

1 研究區概況及數據源

1.1 研究區概況

研究區位于廣西防城港市的北侖河口國家級自然保護區,地理坐標范圍為 E108°02'~108°16',N21°28'~21°37',是一個以紅樹林生態系統為保護對象的自然保護區。該保護區南瀕北部灣,西端與越南交界,自西向東跨越北侖河口、萬尾島和珍珠灣,海岸線全長105 km;有河口海岸、開闊海岸等地貌類型,屬南亞熱帶海洋性季風氣候區。

1.2 數據源

本文研究使用的數據包括Landsat5 TM圖像數據,成像時的潮位資料、氣溫資料(表1),以及紅樹林樣點的野外調查數據。TM圖像通過中國科學院計算機網絡信息中心的國際科學數據鏡像網站(http://datamirror.csdb.cn)下載獲得,成像日期分別為2006年9月12日和2006年10月30日,軌道號/行號為125/45。

表1 TM圖像數據及成像時的潮位和氣溫資料Tab.1 Data list of TM images and the data of tide and air temperature at the time of the image acquisition

2007年11月上旬,作者對研究區內的紅樹林進行了實地野外調查,并使用GPS記錄了各樣點的地理坐標。由于該地區的光學遙感圖像難以獲取,導致本文中圖像獲取與實地野外調查的時間相差1 a。盡管遙感圖像獲取與實地野外調查不同步,但紅樹林為常綠闊葉林,其冠層隨季節和年際變化很小,因此本文野外調查資料可以滿足驗證紅樹林遙感信息提取結果的要求。

2 研究方法

2.1 數據預處理

數據預處理包括對TM圖像的輻射校正和幾何糾正。原始圖像的像元灰度值為0~255的能量相對值,本文對TM圖像進行了輻射定標和大氣校正,Landsat5的輻射定標系數采用了Chander等[18]的研究結果,大氣校正則通過6S程序[19]來完成,輻射校正后獲得各像元的地表反射率(由于本文對熱紅外波段的輻射能量值只考慮其相對大小,因此未對TM6進行輻射校正)。

幾何糾正以1∶5萬地形圖為參考數據,通過選取大型橋梁與河流交叉點、道路交叉口等固定目標作為地面控制點來完成,配準誤差小于0.5個像元。將TM圖像中的熱紅外波段數據的空間分辨率重采樣為30 m,對所有的圖像均進行了UTM投影糾正。

2.2 溫濕度指數

鑒于紅樹林的濕地生長環境,本文基于TM圖像提出了能充分反映紅樹林特點的溫濕度指數(TMI)。該指數能有效地綜合地物的溫度及濕度信息,其定義為

式中SWIR和TIR分別為短波紅外波段反射率和熱紅外波段亮溫。對于TM圖像數據,SWIR=TM5+TM7,主要反映地物的濕度特征;TIR=TM6,主要反映地物的溫度特征。

2.3 分類特征提取與分類算法

分類特征的選擇和提取是地物識別和分類的關鍵。本文采用歸一化均值距離[20]來選擇和評價分類特征。歸一化均值距離通過2個類別的均值距離和標準差來定量評價它們的可分性,距離越大,可分性越好,同時也表明該分類特征較好。其計算公式為

式中:d為歸一化均值距離;μ1和μ2分別為某分類特征的2類樣本區域的均值;σ1和σ2分別為某分類特征的2類樣本區域的標準差。

本文采用光譜角度制圖(spectral angle mapping,SAM)監督分類方法提取紅樹林空間分布信息,將研究區地物分為紅樹林和非紅樹林,并基于ENVI軟件中IDL編程語言實現了SAM監督分類。

3 結果與分析

3.1 典型地物光譜特征

圖1 研究區典型地物反射光譜特征曲線(圖(b)為圖(a)中450~700 nm譜段的放大圖)Fig.1 Reflectance spectral curves of typical objects in the study area

圖1為研究區內紅樹林、森林、農作物(主要為水稻和甘蔗)、城鎮、裸土、灘涂和水體等7種典型地物的反射光譜曲線。其中,植被與非植被之間的光譜差異較大,例如處于短波紅外波段的TM5和TM7可以明顯地將非植被(城鎮、裸土、水體等)與植被區分開來;而不同植被類型之間的光譜極其相似,存在較大的光譜重疊區域;但是,紅樹林與森林、農作物等非紅樹林之間仍然存在一定的光譜差異。從光譜帶重疊的情況來看,紅樹林光譜在可見光和近紅外波段與其他植被類型存在嚴重的重疊,而在對水分含量高度敏感的短波紅外波段則無明顯重疊。因此,短波紅外波段是區分紅樹林與陸地植被的有效波段。

進一步考察水體—植被的混合像元的反射光譜特征可以發現,水體—植被混合像元的反射光譜曲線與紅樹林的極其相似,尤其在短波紅外波段幾乎重疊(圖2)。因此,水體—植被混合像元是造成紅樹林識別精度難以提高的主要原因,僅僅基于反射波段的信息難以對紅樹林與水體—植被混合像元進行有效區分。

圖2 紅樹林與水體—植被混合像元反射光譜特征曲線Fig.2 Spectrum curves of mangrove and water-vegetation mixed pixels

3.2 典型地物溫度特征

為了考察紅樹林、森林、農作物、城鎮、裸土、海水、內陸水體、沿海養殖水體—植被混合像元、內陸水體—植被混合像元等9種典型地物的溫度特征,采用目視判讀方法分別對其進行隨機取樣。由于TM6波段的DN值與像元亮溫之間的函數關系為單調遞增,而本文只考慮像元亮溫的相對大小,因此直接根據TM6波段的DN值進行統計分析(為表達方便,下文均直接借用“亮溫”來表述),統計結果如表2所示。

從表2可以看出,2006年9月12日各地物類型亮溫由小到大依次為:海水、紅樹林、沿海養殖水體—植被混合像元、農作物、內陸水體、內陸水體—植被混合像元、森林、城鎮和裸土。可以發現,沿海海域亮溫普遍低于內陸區域亮溫;且無論在海域還是在內陸,水體均低于同一區域中的植被以及植被—水體的混合像元,這是由于水的比熱大、熱慣量高、白天升溫較慢所致。因此,從海邊向內陸氣溫逐漸升高[21];處于沿海區域的紅樹林亮溫明顯低于森林、內陸水體—植被混合像元的亮溫,但與沿海養殖水體—植被混合像元的亮溫幾乎無差異。

2006年10月30日各地物類型亮溫由小到大依次為:海水、紅樹林、內陸水體、沿海養殖水體—植被混合像元、內陸水體—植被混合像元、森林、農作物、城鎮和裸土。與2006年9月12日的各地物類型亮溫相比,除紅樹林和農作物外,其他地物之間亮溫相對差異無明顯變化。2006年10月30日紅樹林處于高潮位(表1),紅樹林幾乎完全浸泡于海水中,因而紅樹林的亮溫與海水幾乎無差異。農田區域亮溫的大幅提高是因為10月30日已經進入收割季節,大部分耕地變為冬閑地。另外,2006年9月上旬后期至中旬前期,廣西境內出現寒露風天氣,并出現了較明顯的降溫過程(表1),使9月12日各地物的亮溫明顯低于10月30日。

亮溫雖然忽略了各地物類型之間比輻射率的差異,但仍能反映各地物類型之間的溫度特征差異。因此,利用溫度特征可以充分提高紅樹林與森林、農作物的可分性;在結合紅樹林區域潮位資料的基礎上,利用溫度特征還可以進一步提高紅樹林與水體和植被混合像元的可分性。

表2 典型地物溫度特征統計Tab.2 Statistics of typical object temperature characteristics

3.3 分類特征選擇與提取

紅樹林的識別精度主要取決于與陸地植被的區分度,尤其是與水體—植被混合像元的可分性。提取能有效提高紅樹林與陸地植被及水體—植被混合像元的可分性的分類特征(或指數)是解決紅樹林識別精度問題的關鍵。

本文基于TM圖像數據,采用歸一化均值距離作為指標來定量評價紅樹林與其他地物之間的可分性(表3)。在區分紅樹林與其他植被時,對于單波段數據而言,短波紅外波段的可分性最好;而熱紅外波段TM6則可以較好地將紅樹林與水體—植被混合像元分開。

表3 紅樹林與其他地物的區分度Tab.3 Separability between mangrove and other object

紅樹林的特殊生長環境不僅使紅樹林的光譜具有典型植被的特征,同時還反映了影響其生長的濕地背景信息,這使紅樹林與水體、城鎮、裸土等非植被地物非常容易區分;但另一方面也同時導致紅樹林與陰影區陸地植被及植被—水體混合像元的可分性較差。在以往的紅樹林遙感識別研究中,溫度信息往往不受重視,在分析時一般忽略了熱紅外波段數據[3,7]。本文分析認為,雖然熱紅外波段數據的空間分辨率較低,但結合高潮位圖像的溫度特征能夠明顯地提高紅樹林與水體—植被混合像元的區分度(詳見表3中的 TM6H)。另外,通過對原始的反射波段信息和熱紅外波段信息的進一步整合優化發現,本文提出的溫濕度指數((TM5+TM7)TM6)能明顯地提高紅樹林與森林、農作物、裸土和城鎮的可分性,其原因是溫濕度指數揭示了濱海濕地背景與陸地背景的差別。對比已有文獻中通常選取的分類特征發現,(TM5 - TM7)/(TM5+TM7)[4]和 TM4/TM5(注:文獻[6]采用的是TM5/TM4)只能較好地對紅樹林與城鎮、裸土進行區分,NDVI[17]只能較好對紅樹林與城鎮、水體進行區分。

基于上述分析,最終選擇2組分類特征進行紅樹林識別研究:第一組為TM4L(下標L和H分別表示低潮位和高潮位,在本文中分別對應9月12日和10 月 30 日),TM6L,(TM5L+TM7L)TM6L,TM6H,(TM5H+TM7H)TM6H和TM4L/TM5L;第二組作為參照組,用來與本文識別結果進行比較分析,選取的分類特征為其他文獻中常采用的 TM3L,TM4L,TM5L,TM7L,TM3L/TM5L,TM4L/TM5L和 NDVIL。

3.4 紅樹林遙感識別結果分析

首先采用上述兩組分類特征,將地物分成紅樹林和非紅樹林2類,并基于目視判讀方法從TM圖像中提取出441個典型的紅樹林訓練樣本,取其均值波譜作為參考波譜;然后采用SAM監督分類方法進行分類,并對計算后所得訓練樣本的光譜夾角反復進行分析,分別確定2組分類特征的角度閾值(本文用角度表示光譜夾角),第一組和第二組分類特征的角度閾值分別為12°和14°。

以地面實測數據為基礎,采用目視判讀方法從TM圖像中提取519個紅樹林樣本和3880個非紅樹林樣本(包括水體、陸地植被、水體—植被混合像元、城鎮、裸土),對上述2組分類的紅樹林識別結果進行驗證。紅樹林識別精度評價結果詳見表4,兩組分類的紅樹林識別結果詳見圖3。

表4 紅樹林分類結果評價Tab.4 Classification assessment of mangrove

圖3 2006年廣西壯族自治區防城港市北侖河口國家級自然保護區紅樹林識別結果Fig.3 Mangrove map of Beilunhekou National Nature Reserve Area,Fangchenggang City,Guangxi Zhuang Autonomous Region in 2006

由于第一組分類(圖3(a))中將紅樹林的生境特點和TM圖像特點進行了有機結合,充分考慮了像元濕度和溫度特征,選擇和提取了TM6和TMI等能將紅樹林與其他植被類型(陸地植被、水體—植被混合像元)有效區分的分類特征,與未采用熱紅外波段信息的第二組分類(圖3(b))相比,其分類精度有了明顯的提高(Kappa系數提高了0.14,錯分率降低了19.9%)。通過對比分析發現,兩組分類的精度差異主要體現在將植被—水體(包括沿海養殖水體、內陸水體)混合像元錯分為紅樹林(圖3(b))方面,第一組分類中采用了可分性較高的TM6和TMI分類特征,因而有效地降低了紅樹林的錯分率。

4 結論

1)在對水分含量高度敏感的短波紅外波段,紅樹林像元的反射光譜特征與陸地植被(如森林、農作物等)存在較大的差異,而紅樹林與水體—植被混合像元的反射光譜特征極其相似。因此,短波紅外波段是區分紅樹林與陸地植被的有效波段,而僅僅基于反射波段的信息難以滿足紅樹林遙感信息提取的要求。

2)在紅樹林的遙感識別中,通過充分利用紅樹林的特殊生境特點,并基于反射波段特征和熱紅外波段特征提出了溫濕度指數(TMI)。研究結果表明,溫濕度指數能有效地提高紅樹林的識別精度;進一步結合不同潮位信息,則能明顯提高紅樹林與水體—植被混合像元的識別精度。

3)雖然引入熱紅外波段信息可以降低紅樹林與植被—水體混合像元的錯分率,有效提高紅樹林的識別精度,但是仍有部分紅樹林被錯分,這主要是因為TM圖像的空間分辨率偏低(尤其是熱紅外波段的空間分辨率只有120 m)所致。因此,嘗試利用更高空間分辨率的光學圖像(如SPOT4,SPOT5圖像)和熱紅外圖像(如Landsat7 ETM+)來識別紅樹林是下一步的研究重點,預計其結果將會大幅度減少水體—植被混合像元數,進一步提高紅樹林的識別精度。

志謝:中國科學院計算機網絡信息中心國際科學數據鏡像網站 (http://datamirror.csdb.cn)為本文研究免費提供了TM數據,在此表示衷心感謝。

[1]Blasco F,Saemger P,Janodet E.Mangroves as Indicators of Coastal Change[J].Catena,1996,27(3/4):167 -178.

[2]Liu K,Li X,Shi X,et al.Monitoring Mangrove Forest Changes Using Remote Sensing and GIS Data with Decision-tree Learning[J].Wetlands,2008,28(2):336 -346.

[3]Giri C,Pengra B,Zhu Z,et al.Monitoring Mangrove Forest Dynamics of the Sundarbans in Bangladesh and India Using Multi-temporal Satellite Data from 1973 to 2000 [J].Estuarine,Coastal and Shelf Science,2007,73(1/2):91 -100.

[4]Chaudhury M U.Digital Analysis of Remote Sensing Data for Monitoring the Ecological Status of the Mangrove Forests of Sunderbans in Bangladesh[C]//Proceedings of the 23rd International Symposium on Remote Sensing of the Environment,1990,1:493 - 497.

[5]Aschbacher J,Ofren R S,Delsol J P,et al.An Integrated Comparative Approach to Mangrove Vegetation Mapping Using Advanced Remote Sensing and GIS Technologies:Preliminary Results[J].Hydrobiologia,1995,295(1/3):285 -294.

[6]Green E P,Clear C D,Mum P J,et al.Remote Sensing Techniques for Mangrove Mapping[J].International Journal of Remote Sensing,1998,19(5):935 -956.

[7]Andriamparany R,Francois F.Dynamics of Mangrove Forests in the Mangoky River Delta,Madagascar,Under the Influence of Natural and Human Factors[J].Forest Ecology and Management,2010,259(6):1161-1169.

[8]Blasco F,Aizpuru M,Gers C.Depletion of the Mangroves of Continental Asia[J].Wetlands Ecology and Management,2001,9(3):245-256.

[9]Thu P M,Pppulus J.Status and Changes of Mangrove Forest in Mekong Delta:Case Study in Tra Vinh,Vietanam[J].Estuarine Coastal and Shelf Science,2007,71(1/2):98 -109.

[10]Conchedda G,Laurent D,Philippe M.An Object-based Method for Mapping and Change Analysis in Mangrove Ecosystems[J].ISPRS Journal of Photogrammetry and Remote Sensing,2008,63(5):578-589.

[11]Vaiphasa C,Andrew S K,Willem F B.A Post-classifier for Mangrove Mapping Using Ecological Data[J].ISPRS Journal of Photogrammetry and Remote Sensing,2006,61(1):1 -10.

[12]肖海燕,曾 輝,昝啟杰,等.基于高光譜數據和專家決策法提取紅樹林群落類型信息[J].遙感學報,2007,11(4):531-537.Xiao H Y,Zeng H,Zan Q J,et al.Decision Tree Model in Extraction of Mangrove Community Information Using Hyperspectral Image Data[J].Journal of Remote Sensing,2007,11(4):531 -537(in Chinese with English Abstract).

[13]Wang L,Sousa W P,Gong P,et al.Comparison of IKONOS and QuickBird Images for Mapping Mangrove Species on the Caribbean Coast of Panama [J].Remote Sensing of Environment,2004,91(3/4):432-40.

[14]Kovacsa J M,Wang J F,Flores- Verdugoc F.Mapping Mangrove Leaf Area Index at the Species Level Using IKONOS and LAI-2000 Sensors for the Agua Brava Lagoon Mexican Pacific[J].Estuarine,Coastal and Shelf Science,2005,62(1/2):377 -384.

[15]Proisy,Coutero P,Fromard F.Predicting and Mapping Mangrove Biomass from Canopy Grain Analysis Using Fourier-based Textural Ordination of IKONOS Images[J].Remote Sensing of Environment,2007,109(3):379 -392.

[16]Saleh A.Assessment of Mangrove Vegetation on Abu Minqar Island of the Red Sea[J].Journal of Arid Environments,2007,68(2):331-336.

[17]Jensen J R,Ramset E,Davis B A,et al.The Measurement of Mangrove Characteristics in South-west Florida Using SPOT Multispectral Data[J].Geocartography International,1991,6(2):13 -21.

[18]Chander G,Markham B.Revised Landsat-5 TM Radiometric Calibration Procedures and Postcalibration Dynamic Ranges[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(11):2674-2677.

[19]Vetmote E F,Tanre D,Deuze J L,et al.Second Simulation of the Satellite Signal in the Solar Spectrum,6s:An Overview[J].IEEE Transactions on Geoscience and Remote Sensing,1997,35(3):675-686.

[20]Swain P H,Davis S M.Remote Sensing:The Quantitative Approach[M].New York:McGrowHill Inc,1978.

[21]朱萊茵,許映軍,崔維佳,等.渤海灣西部沿岸地區氣溫特征的觀測研究[J].氣象科學,2009,29(5):694 -699.Zhu L Y,Xu Y J,Cu W J,et al.The Observation Research of Temperature Impacted by the Sea Land Breeze in the Coastal Area of Western Bohai Bay [J].Scientia Meteorologica Sinica,2009,29(5):694-699(in Chinese with English Abstract).

Application of the Temperature-Moisture Index to the Improvement of Remote Sensing Identification Accuracy of Mangrove

ZHANG Xue-hong1,TIAN Qing-jiu2
(1.School of Remote Sensing,Nanjing University of Information Science & Technology,Nanjing 210044,China;2.International Institute for Earth System Science,Nanjing University,Nanjing 210093,China)

The identification accuracy of mangrove by using TM reflective bands is always low due to the similarity of spectra between mangrove and land vegetation,especially water- vegetation mixed pixels.Based on reflective and thermal infrared information in the TM images of different tide levels,the authors proposed temperature -moisture index(TMI).The analytical results show that the thermal infrared band and TMI can obviously improve the separability between mangrove and other objects based on the tide level information.The thermal infrared band and TMI can also significantly increase the classification accuracy of mangrove by using spectral angle mapping(SAM) supervised classification method in comparison with the classification features employed by other researchers.The Kappa coefficient increases by 0.14 ,and the commission error of mangrove class decreases by 19.9%,suggesting that the remote sensing identification accuracy of mangrove can be improved by using the information of tide level,thermal infrared band and TMI.

mangrove;tide level;TM;thermal infrared;temperature-moisture index(TMI)

TP 751.1;TP 753

A

1001-070X(2012)03-0065-06

2011-11-24;

2012-01-02

國家自然科學基金項目(編號:41201461)、國防科技工業民用科研技術研究項目(編號:2006A100602)和江蘇高校優勢學科建設工程資助項目共同資助。

10.6046/gtzyyg.2012.03.13

張雪紅(1980-),男,博士,講師,主要從事植被生態遙感方面的研究。E-mail:zxhbnu@126.com。

(責任編輯:劉心季)

猜你喜歡
紅樹林分類特征
走過紅樹林
歌海(2024年6期)2024-03-18 00:00:00
藏著寶藏的紅樹林
分類算一算
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
分類討論求坐標
神奇的紅樹林
數據分析中的分類討論
走過紅樹林
歌海(2018年4期)2018-05-14 12:46:15
教你一招:數的分類
主站蜘蛛池模板: 99999久久久久久亚洲| 国产成人精品免费视频大全五级| 欧美成人精品一级在线观看| 免费中文字幕一级毛片| 在线观看91香蕉国产免费| 欧美在线一二区| 日韩精品一区二区三区中文无码| 成年人久久黄色网站| 国产尹人香蕉综合在线电影| 黄色网在线| 无码区日韩专区免费系列| 中文字幕啪啪| 亚洲激情99| 午夜无码一区二区三区| 自拍欧美亚洲| 日韩av无码精品专区| 日韩欧美综合在线制服| 国内黄色精品| 国产幂在线无码精品| 免费毛片在线| 亚洲美女视频一区| 曰AV在线无码| 亚洲乱码精品久久久久..| 久久96热在精品国产高清| 亚洲精品国产首次亮相| 午夜a级毛片| 狠狠色丁香婷婷| 国产午夜小视频| 在线观看国产小视频| 波多野结衣在线se| 亚洲三级电影在线播放| 亚洲精品动漫| 亚洲国产看片基地久久1024| 91成人精品视频| 特级毛片8级毛片免费观看| 高清精品美女在线播放| 乱系列中文字幕在线视频| 国产精品美女免费视频大全| 国产成人亚洲毛片| 久久精品中文字幕免费| 久久午夜夜伦鲁鲁片不卡| 久久综合伊人77777| 国产91视频观看| 国产成人做受免费视频| 99视频国产精品| 久久精品国产一区二区小说| 狠狠躁天天躁夜夜躁婷婷| 国产小视频在线高清播放 | 中文字幕在线看| 日本国产在线| 天堂av高清一区二区三区| 成人精品免费视频| 国产成人久久综合一区| 一级片免费网站| 欧美性精品| 直接黄91麻豆网站| 久久久久亚洲av成人网人人软件 | 久无码久无码av无码| 国内老司机精品视频在线播出| 亚洲第一视频区| 亚洲AV电影不卡在线观看| 欧美日韩国产系列在线观看| 亚洲无线一二三四区男男| 久久网综合| 国产精选小视频在线观看| 无码网站免费观看| 四虎AV麻豆| 国产精品99久久久久久董美香| P尤物久久99国产综合精品| 国产精品亚洲va在线观看| 8090午夜无码专区| 精品一区二区三区自慰喷水| 国产91丝袜在线播放动漫 | www.精品国产| 五月激激激综合网色播免费| 久99久热只有精品国产15| 巨熟乳波霸若妻中文观看免费| 美女一区二区在线观看| 久久黄色毛片| 亚洲成年人网| 毛片基地视频| 中文字幕在线不卡视频|