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

山西省環(huán)境星數(shù)據(jù)纓帽變換參量的提取研究

2016-12-15 07:28:45崔鵬艷李大成
關(guān)鍵詞:研究

崔鵬艷,李大成,馬 綽

(太原理工大學(xué) 礦業(yè)工程學(xué)院,太原 030024)

?

山西省環(huán)境星數(shù)據(jù)纓帽變換參量的提取研究

崔鵬艷,李大成,馬 綽

(太原理工大學(xué) 礦業(yè)工程學(xué)院,太原 030024)

針對(duì)環(huán)境一號(hào)衛(wèi)星(HJ-1)CCD數(shù)據(jù)纓帽變換參數(shù)的提取問題,提出了傳統(tǒng)的纓帽參數(shù)提取方法與光譜分析技術(shù)、施密特正交化算法相結(jié)合的HJ-1 CCD數(shù)據(jù)的纓帽變換參數(shù)提取算法,利用該方法對(duì)山西省全區(qū)2009-2013年間處于植被生長(zhǎng)季節(jié)的HJ-1 CCD數(shù)據(jù)進(jìn)行實(shí)驗(yàn)驗(yàn)證。結(jié)果表明,利用提取得到的纓帽參數(shù)對(duì)研究區(qū)的上述數(shù)據(jù)進(jìn)行纓帽變換處理后,圖像特征得到顯著增強(qiáng),亮度分量、綠度分量以及濕度分量對(duì)于地表特征的表達(dá)更為合理與清晰。該方法為區(qū)域條件下地表特征的識(shí)別與表達(dá)提供有效的技術(shù)途徑。

光譜統(tǒng)計(jì)分析;施密特正交化;HJ-1 CCD;纓帽參數(shù);環(huán)境星數(shù)據(jù);山西省

我國(guó)自行研制的環(huán)境與災(zāi)害監(jiān)測(cè)預(yù)報(bào)系統(tǒng)一號(hào)衛(wèi)星的A、B星座(以下簡(jiǎn)稱HJ-1 A/B)于2008年9月6日在太原衛(wèi)星發(fā)射中心成功發(fā)射,至今正常運(yùn)行,并在遙感數(shù)據(jù)的大氣校正[1-2]、干旱監(jiān)測(cè)[3]、云檢測(cè)[4]以及環(huán)境質(zhì)量評(píng)價(jià)[5]等應(yīng)用領(lǐng)域發(fā)揮了重要作用。

HJ-1 A/B衛(wèi)星均裝載了相同類型的CCD相機(jī),通過在相同軌道上星下點(diǎn)的對(duì)稱放置以達(dá)到平分視場(chǎng)、并行觀測(cè)的目的;同時(shí)可使組網(wǎng)后CCD相機(jī)的重訪周期僅為2 d,極大地提高了該衛(wèi)星的應(yīng)用潛力。另一方面,A、B星所搭載的CCD相機(jī)只包含4個(gè)常見光譜通道(如表1所示),顯然無(wú)法充分滿足多數(shù)應(yīng)用對(duì)中分辨率遙感數(shù)據(jù)波段數(shù)量的要求。因此,如何針對(duì)HJ-1 A/B CCD相機(jī)的數(shù)據(jù)條件,開發(fā)和設(shè)計(jì)出有效的光譜變換方法,就成為該問題的重要解決途徑。基于多光譜數(shù)據(jù)的光譜變換可以獲得比原有圖像更加結(jié)構(gòu)化的光譜分布特征。常見的光譜變換方法主要有主成分分析(principal components analysis,PCA)和纓帽變換(tasseled cap,TC)。PCA 方法是通過降維,盡可能多地保留原始特征集合中的特征信息量,得到線性無(wú)關(guān)的綜合向量[6];TC變換是一種針對(duì)地表信息提取的影像增強(qiáng)方法,是根據(jù)多光譜遙感中土壤、植被等信息在多維光譜空間中信息分布結(jié)構(gòu)對(duì)圖像做的經(jīng)驗(yàn)性正交變換[7],相較于PCA,它具有更為明確的地學(xué)意義。

TC變換是KAUTH和THOMAS在研究Landsat MSS圖像反映農(nóng)作物和植被生長(zhǎng)過程時(shí)提出的[8],他們對(duì)Landsat MSS進(jìn)行TC參數(shù)的提取和研究時(shí),發(fā)現(xiàn)提取的前3個(gè)分量分別為亮度、綠度和黃度分量。CHRIST和CICONE利用TC變換對(duì)Landsat TM進(jìn)行纓帽參數(shù)提取研究時(shí),發(fā)現(xiàn)其第3個(gè)分量并非表示黃度值,而是濕度值[9]。黃成全等人利用Landsat 7 ETM+數(shù)據(jù)進(jìn)行TC變換的研究時(shí)發(fā)現(xiàn),利用DN值提取的纓帽參數(shù)是不準(zhǔn)確的;而基于表觀反射率提取的亮度、綠度、濕度分量是相互獨(dú)立,并且不受云污染的影響[10]。SCHAAF利用基于Modis 數(shù)據(jù)的纓帽變換來研究植被的周期變換規(guī)律[11]; YARBROUGH發(fā)展了基于Quickbird數(shù)據(jù)的纓帽參數(shù)提取方法研究[12]。雖然當(dāng)前圍繞不同傳感器已開展相應(yīng)的纓帽參數(shù)提取研究,但針對(duì)HJ-1 CCD相機(jī)的纓帽參數(shù)提取仍缺少區(qū)域級(jí)高時(shí)序數(shù)據(jù)的支撐及統(tǒng)計(jì)分析。因此,筆者針對(duì)HJ-1 A/B數(shù)據(jù)的纓帽參數(shù)的提取展開研究,為區(qū)域條件下地表特征的識(shí)別與表達(dá)提供有效的技術(shù)途徑。并利用山西省2009—2013年間處于植被生長(zhǎng)季節(jié)的HJ-1 CCD數(shù)據(jù),為基于其HJ-1 CCD數(shù)據(jù)的纓帽變換處理提供更加精確的轉(zhuǎn)換參數(shù)。

1 研究區(qū)與數(shù)據(jù)預(yù)處理

本研究以山西省全區(qū)為實(shí)驗(yàn)研究區(qū)。研究區(qū)地處我國(guó)黃河中游東岸,屬華北平原西部的黃土高原地帶,地理位置為:北緯34°34′-40°44′,東經(jīng)110°14′-114°33′,總面積15.67×104km2。研究區(qū)與鄰省的自然境界分明,地勢(shì)東北高西南低;區(qū)內(nèi)部起伏不平,河谷縱橫,地貌類型復(fù)雜多樣,有山地、平原、丘陵、臺(tái)地,植被生長(zhǎng)季節(jié)為每年的5月至9月。

本實(shí)驗(yàn)選取覆蓋研究區(qū)2009-2013年間共5 a的HJ-1 A/B CCD數(shù)據(jù)為基礎(chǔ)數(shù)據(jù)。選取數(shù)據(jù)的基本原則為成像質(zhì)量較高且云量覆蓋小于10%,并在上述年份的所有植被生長(zhǎng)季節(jié),逐月?lián)袢「采w研究區(qū)北部和南部且同月獲取時(shí)間接近、鄰月相差超過16 d的影像數(shù)據(jù)。最終篩選得到47幅研究數(shù)據(jù),2009-2013年間的植被生長(zhǎng)季相內(nèi)圖像選取數(shù)量的分布關(guān)系如圖1所示。

圖1 研究區(qū)環(huán)境星數(shù)據(jù)的季相分布Fig.1 Seasonal distribution of the HJ-1 CCD data of the research area

在實(shí)驗(yàn)研究前,需對(duì)所選取的環(huán)境星影像進(jìn)行必要的預(yù)處理。首先,利用ENVI軟件通過輻射定標(biāo)原理,將圖像中原始的無(wú)量綱DN值(Digital Number,DN)數(shù)據(jù)轉(zhuǎn)換為相應(yīng)的輻射亮度值;然后根據(jù)相關(guān)參數(shù)在ENVI中將所得的輻射亮度值進(jìn)一步轉(zhuǎn)換為大氣層頂表觀反射率(Top of Atmosphere,TOA);最后,為了研究的方便,將處理后的TOA值域拉伸至[0,10 000]。

2 基于環(huán)境星數(shù)據(jù)的纓帽系數(shù)提取原理

原始TC變換的模型如式(1)所示,是KAUTH和THOMAS通過對(duì)Landsat MSS圖像的研究提出的。HJ-1 A/B數(shù)據(jù)與Landsat MSS數(shù)據(jù)具有相同的波段數(shù),其對(duì)應(yīng)波段的基本參數(shù)信息如表1所示。分析可知,兩者的波段通道信息相差不大,因而基于Landsat MSS數(shù)據(jù)的TC變換同樣適用于基于HJ-1 A/B數(shù)據(jù)的TC變換研究,且HJ-1 A/B數(shù)據(jù)比Landsat MSS數(shù)據(jù)的空間分辨率高,因此針對(duì)HJ-1 A/B數(shù)據(jù)的應(yīng)用價(jià)值更高。

表1 HJ-1 A/B 數(shù)據(jù)與Landsat MSS數(shù)據(jù)的基本參數(shù)

原始TC變換模型

(1)

式中:X為L(zhǎng)andsat MSS各波段的數(shù)據(jù)信息;U為纓帽變換后提取的數(shù)據(jù)信息;RT為轉(zhuǎn)換矩陣[6]。同理,易得基于環(huán)境星CCD圖像的纓帽變換模型為:

(2)

式中:Xi為變化前第i波段上的TOA值;RT是由纓帽參數(shù)Bi,Gi,Wi構(gòu)成的系數(shù)矩陣(i=1,2,3,4);U為變化后由亮度、綠度及濕度分量組成的矩陣。

式(2)中的亮度分量反應(yīng)了土壤反射率變化的信息,其系數(shù)項(xiàng)(Bi)可以利用干土與濕土的反射光譜差異得到,計(jì)算過程可表達(dá)如下:

(3)

(4)

(5)

式(2)中的綠度分量系數(shù)(Gi)可以提取諸如地面植被等圖像綠度信息,結(jié)合植被與干土的反射光譜差異以及(5)式獲取的亮度分量系數(shù)(Bi),利用式(6)、式(7)求得與亮度分量正交的綠度分量;并根據(jù)式(8)、式(9)將正交化結(jié)果權(quán)值化,從而得到歸一化的綠度分量系數(shù)。計(jì)算過程可表達(dá)如下:

(6)

(7)

(8)

(9)

式(2)中的濕度分量系數(shù)(Wi)可以提取水體等濕度信息,根據(jù)水體與干土的反射光譜差異以及式(5)、式(9)得到的Bi、Gi,利用式(10)、(11)、(12)求得與亮度分量、綠度分量正交的濕度分量系數(shù),并根據(jù)式(13)、(14)將正交結(jié)果權(quán)值化,得到歸一化的濕度分量系數(shù)。其計(jì)算過程如下:

(10)

(11)

(12)

(13)

(14)

亮度分量B、綠度分量G、濕度分量W可以利用提取的亮度分量系數(shù)、綠度分量系數(shù)以及濕度分量系數(shù)與原始TOA值的乘積得到,表達(dá)式如下:

(15)

(16)

(17)

為得到基于實(shí)驗(yàn)研究區(qū)的TC變換的系數(shù)矩陣,需在數(shù)據(jù)預(yù)處理后的TOA圖像上人工提取出干土、濕土、植被以及水體等4類典型地物的光譜統(tǒng)計(jì)信息,本研究通過逐波段統(tǒng)計(jì)地物特征的光譜均值來實(shí)現(xiàn)。

3 實(shí)驗(yàn)結(jié)果與分析

圖2中的a-d分別描述了研究區(qū)內(nèi)的不同地物類別分別在藍(lán)光、綠光、紅光以及近紅外波段上的反射光譜信息隨時(shí)序的變化特征。其中,橫坐標(biāo)為數(shù)據(jù)的采集時(shí)間(季相),縱坐標(biāo)為各種地物特征在相應(yīng)波段上的反射光譜均值;平滑實(shí)線、平滑虛線、三角實(shí)線以及方塊實(shí)線分別代表干土、濕土、植被以及水體類型的地物反射光譜變化曲線。從圖2可見,以上4種地物類別在藍(lán)光波段的反射光譜的強(qiáng)弱關(guān)系為ρdrysoi,ρvegetation>ρwetsoil>ρwater。植被和水體在綠光波段的反射光譜相近,分布范圍為[0,1 000];而濕土的反射光譜分布于[1 000,2 500],干土的反射光譜整體分布于[2 000,4 000],且ρdrysoil>ρwetsoil>ρvegetation,ρwater。植被與水體在紅光波段的反射光譜相近(分布區(qū)間為[400,1 000]),濕土在紅光波段的分布區(qū)間為[1 100,1 800],干土在紅光波段的反射特征最強(qiáng)(分布于[1 300,3 000]);但其波動(dòng)性較大,反射光譜的強(qiáng)弱關(guān)系與綠光波段類似。各類型地物在近紅外波段反射光譜的可區(qū)分性稍差,整體上有:ρdrysoil>ρwetsoil>ρvegetation,ρwater。

圖2 不同地類在各波段上的反射率隨時(shí)序的變化特征Fig.2 Change characteristics of the reflectance of different surface features on various wavelengths along with the time sequence

由于纓帽變換的精度受季相性和區(qū)域性的影響較大[13],筆者所選的實(shí)驗(yàn)研究區(qū)在地理分布上特點(diǎn)是東西向較窄,南北向跨度較大;因此對(duì)南北部的典型地物類別(植被,水體,干土以及濕土)在各波段上的反射光譜(表觀反射率)特征進(jìn)行了統(tǒng)計(jì)分析,了解研究區(qū)南部與北部典型地物特征的反射光譜的一致性程度。

圖3所示為4種典型地物在近紅外波段的南北反射光譜信息,總體上4類地物反射光譜信息在南部與北部分布一致。其中,植被南北數(shù)據(jù)分布一致;水體在2009年6月差異較大;干土于2009年8月差異較大;濕土南北數(shù)據(jù)波段較大,但兩者分布吻合。因此,植被、水體、干土、濕土在研究區(qū)的分布受南北地域影響不大。利用TC參量的提取策略與計(jì)算方法所得到的基于實(shí)驗(yàn)研究區(qū)的纓帽變換系數(shù)的統(tǒng)計(jì)信息如表2所示。

表2 山西省全區(qū)的纓帽變換參量的統(tǒng)計(jì)信息

為了檢測(cè)提取的亮度、綠度以及濕度系數(shù)的應(yīng)用效果,利用提取的纓帽變換系數(shù)獲取其亮度、綠度以及濕度分量。圖4所示是以假彩色合成圖像為基礎(chǔ)進(jìn)行的近紅外波段與亮度分量的對(duì)比分析。其中近紅外圖像上云都呈現(xiàn)亮白色,屬于強(qiáng)反射;水都呈現(xiàn)暗灰色,屬于弱反射。亮度圖像中云的顏色變亮,水的顏色變暗,其中2009年和2012年效果最明顯。從近紅外與亮度分量(圖4右)的定量分析中,虛線為近紅外光譜反射信息,實(shí)線為亮度分量信息;橫坐標(biāo)為拉伸后的地表反射率值或亮度值,縱坐標(biāo)為像元個(gè)數(shù),其統(tǒng)計(jì)結(jié)果分布類似正態(tài)分布,整體上亮度的均值和方差比近紅外的大。其中,2009年和2012年分布最明顯,與定性分析一致。

圖3 研究區(qū)內(nèi)典型地物在南部與北部區(qū)域的TOA反射率對(duì)比Fig.3 The comparison of TOA reflectance between the south and the north of study area about typical surface features

圖4 近紅外圖像與亮度圖像的對(duì)比Fig.4 Comparison between the near-infrared data and the brightness data

歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)應(yīng)用于檢測(cè)植被生長(zhǎng)狀態(tài)、植被覆蓋度和消除部分輻射誤差等[14],可以檢查綠度分量的效果。圖5所示是以假彩色合成圖像為基礎(chǔ)進(jìn)行的近紅外波段與綠度分量、NDVI的對(duì)比分析。近紅外圖像中植被分布呈現(xiàn)暗灰色,在綠度圖像和NDVI圖像中其分布呈亮白色。在NDVI與綠度的散點(diǎn)圖(圖5右)中,2009—2013年的數(shù)據(jù)都是沿?cái)M合直線對(duì)稱分布。其中,2011年、2012年和2013年數(shù)據(jù)分布較集中。

圖6所示是以假彩色合成圖像為基礎(chǔ)進(jìn)行的近紅外與濕度影像的對(duì)比分析。以水為例,從定性角度分析:水在近紅外影像中為黑色,在濕度影像中為白色,水體變得發(fā)亮;從定量角度分析,變換前水體的TOA值都比圖像整體的TOA值小,而變換后水體的濕度值都比圖像整體的濕度值大。

圖5 NDVI圖像與綠度圖像的線性相關(guān)性對(duì)比Fig.5 Comparison between the NDVI data and the greenness data

圖6 近紅外圖像與濕度圖像對(duì)比信息Fig.6 Comparison between the near-infrared data and the wetness data

4 結(jié)束語(yǔ)

本文以覆蓋山西省全區(qū)2009-2013年間處于植被生長(zhǎng)季節(jié)(5-9月)的HJ-1 A/B CCD為數(shù)據(jù)源,利用類別的光譜分析技術(shù)與施密特正交化算法進(jìn)行研究區(qū)纓帽變換參數(shù)的提取,初步得到以下結(jié)論:

1) 纓帽變換受季相性的影響。本實(shí)驗(yàn)利用植被生長(zhǎng)季節(jié)(5-9月)數(shù)據(jù)所提取的纓帽變換系數(shù),適用于獲取于該季相的環(huán)境星數(shù)據(jù)。

2) 纓帽變換受區(qū)域性的影響。研究區(qū)在地理分布上特點(diǎn)是:東西向窄,南北向?qū)?干土、濕土、植被以及水體4類地物特征的反射光譜信息隨區(qū)域差異(南北部)變化很小,因而所提取的纓帽系數(shù)有望適用于整個(gè)地理研究區(qū)。

3) 纓帽系數(shù)對(duì)于亮度、綠度以及濕度分量的生成效果。實(shí)驗(yàn)分析可知,亮度分量對(duì)于圖像亮度特征的描述更加豐富;綠度分量與NDVI圖像一致性較高,且兩者基本呈線性分布;濕度分量對(duì)于水體等“濕度”特征的描述更加直觀,且其灰度特征更加集中。因此,上述纓帽分量對(duì)于地表特征的后處理操作(如分類)能夠提供更好的結(jié)構(gòu)化光譜分布特征及特征識(shí)別環(huán)境。

[1] 劉其悅.基于環(huán)境星數(shù)據(jù)可見光、近紅外波段的大氣校正方法研究[D].焦作:河南理工大學(xué),2010.

[2] 孫源,顧行發(fā),余濤,等.環(huán)境星CCD數(shù)據(jù)大氣校正研究[J].國(guó)土資源遙感,2010(4):6-9.

[3] 鐘仕全,羅永明,莫建飛,等.環(huán)境減災(zāi)衛(wèi)星數(shù)據(jù)在干旱監(jiān)測(cè)中的應(yīng)用[J].中國(guó)農(nóng)業(yè)氣象,2011,32(4):593-597.

[4] LI D,GE Y,WANG B.Automated retrieval of cloud masks from the HJ-1 WVC imagery[J].Selected Topics in Applied Earth Observations and Remote Sensing,2014(9):3732-3741.

[5] 劉瑞,王世新,周藝,等.基于遙感技術(shù)的縣級(jí)區(qū)域環(huán)境質(zhì)量評(píng)價(jià)模型研究[J].中國(guó)環(huán)境科學(xué),2012,32(1):181-186.

[6] 李建林.一種基于PCA的組合特征提取文本分類方法[J].計(jì)算機(jī)應(yīng)用研究,2013,30(8):2398-2401.

[7] 湯國(guó)安,張友順,劉詠梅,等.遙感數(shù)字圖像處理[M].北京:科學(xué)出版社,2004.

[8] KAUTH R J,THOMAS G S.The tasselled cap-a graphic description of the spectral-temporal development of agricultural crops as seen by landsat[C]∥LARS Symposia.Proceedings of the Symposium of Machine Processing of Remotely Sensed Data. Indiana:The Laboratory for Applications of Remote Sensing,1976:41-51.

[9] CRIST E,CICONE R.A physically-based transformation of Thematic Mapper data-The TM tasselled cap[J].IEEE Transactions on Geoscience and Remote Sensing,1984,22(3):256-263.

[10] HUANG C Q,WYLI B,YANG L,et al.Derivation of a tasseled cap transformation based on Landsat 7 at satellite reflectance[J].International Journal of Remote Sensing,2002,23(8):1741-1748.

[11] ZHANG X,SCHAAF C B,FRIEDL M A,et al.MODIS tasseled cap transformation and its utility[J].Geoscience and Remote Sensing Symposium,2002(2):1063-1065.

[12] YARBROUGH L D,EASSON G,KUSZMAUL J S.QuickBird 2 tasseled cap transform coefficients:a comparison of derivation methods[J].Global Priorities in Land Remote Sensing,2005(16):23-27.

[13] EVA I,ALISTAIR L,FILIP L,et al.Orthogonal transformation of segmented SPOT5 images: seasonal and geographical dependence of the tasselled cap parameters[J].Photogrammetric Engineering and Remote Sensing,2008,74(11):1351-1364.

[14] JUSTICE C.Analysis of the phonology of global vegetation using meteorological satellite data[J].International Journal of Remote Sensing,1985(6):1271-1318.

(編輯:龐富祥)

Extraction of Tasselled Cap Transformation Parameters of Environmental Satellite Data of Shanxi

CUI Pengyan,LI Dacheng,MA Chao

(CollegeofMiningEngineering,TaiyuanUniversityofTechnology,Taiyuan030024,China)

On the basis of traditional method of extraction of Tasselled Cap transformation parameters, this research put forward the extraction of Tasselled Cap transformation parameters for HJ-1 A/B CCD data by combining the technology of spectral analysis based on category and the means of Schmidt orthogonalization. This method was used to extract the Tasselled Cap transformation parameters and verify the accuracy of the HJ-1 CCD date during the vegetation growing season from 2009-2013 in Shanxi. The results show that after transforming the study area of HJ-1 A/B CCD data with the Tasselled Cap transformation parameters, the feature of the image has been significantly enhanced, and the luminance component, greenness component and humidity component for expression of the surface characteristics are more reasonable and clearer. The extracted information has no much difference between the south and the north, so its Tasselled Cap transformation parameters can apply to the whole Shanxi Province during the vegetation growing season. It provides an effective way to identify and express the surface features in the regional conditions.

spectral statistical analysis;schmidt orthogonalization;HJ-1 CCD;tasselled cap transformation parameters;environmetal satellite data;Shanxi province

1007-9432(2016)03-0314-07

2015-10-26

科技部863計(jì)劃基金資助項(xiàng)目:全球地表覆蓋遙感制圖與關(guān)鍵技術(shù)研究(2009AA122002)

崔鵬艷(1989-),女,河南洛陽(yáng)人,碩士生,主要從事遙感地質(zhì)研究,(E-mail)cuipengyancehui@163.com

李大成,男,講師,主要從事遙感信息處理及地學(xué)應(yīng)用研究,(E-mail)lidacheng@tyut.edu.cn

TP79

A

10.16355/j.cnki.issn1007-9432tyut.2016.03.008

猜你喜歡
研究
FMS與YBT相關(guān)性的實(shí)證研究
2020年國(guó)內(nèi)翻譯研究述評(píng)
遼代千人邑研究述論
視錯(cuò)覺在平面設(shè)計(jì)中的應(yīng)用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
關(guān)于遼朝“一國(guó)兩制”研究的回顧與思考
EMA伺服控制系統(tǒng)研究
基于聲、光、磁、觸摸多功能控制的研究
電子制作(2018年11期)2018-08-04 03:26:04
新版C-NCAP側(cè)面碰撞假人損傷研究
關(guān)于反傾銷會(huì)計(jì)研究的思考
焊接膜層脫落的攻關(guān)研究
電子制作(2017年23期)2017-02-02 07:17:19
主站蜘蛛池模板: 日韩av手机在线| 国产黑丝视频在线观看| 日本欧美一二三区色视频| 国产91麻豆免费观看| 亚洲人成网站18禁动漫无码| 无码一区二区三区视频在线播放| Aⅴ无码专区在线观看| 91精品视频在线播放| 啪啪永久免费av| 国产日本视频91| 亚洲V日韩V无码一区二区| 中文字幕 91| 人妻中文字幕无码久久一区| 97精品伊人久久大香线蕉| 亚洲精品在线91| www.av男人.com| 爽爽影院十八禁在线观看| 国内老司机精品视频在线播出| 国产99久久亚洲综合精品西瓜tv| 国产区免费| 国产亚洲精品自在久久不卡| 真实国产精品vr专区| 欧美自拍另类欧美综合图区| 国产小视频免费| 精品视频一区二区三区在线播| 日韩欧美91| 国产幂在线无码精品| 国产三级韩国三级理| 中文天堂在线视频| 欧美有码在线| 国产精品无码制服丝袜| 色窝窝免费一区二区三区 | 国产一级精品毛片基地| 在线国产毛片| 久久综合婷婷| AⅤ色综合久久天堂AV色综合| 欧洲极品无码一区二区三区| 国产高清在线观看91精品| 国产正在播放| 国产91丝袜| 亚洲欧美极品| 国内精品久久人妻无码大片高| 亚洲第一av网站| 国产激情第一页| 又爽又大又黄a级毛片在线视频| 亚洲AⅤ永久无码精品毛片| 色偷偷一区二区三区| 国产精品毛片在线直播完整版| 国产精品无码AV片在线观看播放| 波多野结衣一区二区三区四区视频| 欧美日韩国产在线播放| 日韩色图区| 中文字幕乱码中文乱码51精品| 在线观看91精品国产剧情免费| 亚洲国产日韩一区| 最新亚洲av女人的天堂| 岛国精品一区免费视频在线观看 | 国产成人亚洲毛片| 欧美激情伊人| 被公侵犯人妻少妇一区二区三区| 国产香蕉在线视频| 欧美一区二区三区香蕉视| 精品自拍视频在线观看| 亚洲欧洲日产无码AV| 高清欧美性猛交XXXX黑人猛交| 在线观看亚洲国产| 丝袜高跟美脚国产1区| 在线观看国产一区二区三区99| 亚洲国产欧美目韩成人综合| 98精品全国免费观看视频| www.99在线观看| 92午夜福利影院一区二区三区| 国产精品美女网站| 9966国产精品视频| 在线国产91| 国内精品九九久久久精品| 欧美区一区二区三| 国产交换配偶在线视频| 中文字幕欧美日韩高清| 成年A级毛片| 在线国产欧美| 亚洲欧美日韩中文字幕在线|