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

基于DINEOF方法的水色遙感數據的重構研究
——以黃、渤海區域為例

2014-08-01 09:33:16王躍啟劉東艷
遙感信息 2014年5期
關鍵詞:方法

王躍啟,劉東艷

(1.中國科學院煙臺海岸帶研究所 海岸帶環境過程與生態修復重點實驗室,山東 煙臺 264003;2.中國科學院大學,北京 100049)

1 引 言

海洋水色遙感具有成本低、數據量大,可實現大范圍連續觀測等優點,并且與傳統監測方法有很好的互補性,在海洋生態系統的研究中扮演著越來越重要的作用[1-4]。近年來,海洋水色遙感技術發展迅速,尤其是SeaWiFS(Sea-viewing Wide Field-of-view Sensor)和MODIS(Moderate Resolution Imaging Spectrotadiometer) 等傳感器的業務化運行,不僅為海洋研究提供了長時間序列的水色遙感數據產品,而且實現了多傳感器的同步、互補觀測,為海洋生態系統的長期研究提供了豐富的數據源。

隨著海洋水色遙感數據的不斷積累,目前海洋水色遙感數據產品的應用仍面臨著兩個關鍵的問題:一方面,受云覆蓋、大氣校正、傳感器技術問題等的影響,獲取的數據面臨著較高的缺失問題,需要有效的空間插值方法來提高數據的時空完整性;另一方面,不同傳感器的運行時間和壽命不同,需要有效的融合方法將不同的傳感器數據進行準確的融合和承接,最大限度地獲得長時間序列的完整數據,以滿足對海洋環境的長期連續觀測要求[5]。為了解決上述問題,國內外學者發展了一系列的插值和融合方法。傳統的海洋水色遙感數據的插值方法主要包括地統計插值法(geostatistical filling)[6-7]、最優插值法(Optimal Interpolation,OI)[8]、奇異譜分析法 (Singular Spectrum Analysis,SSA)[9]等;而遙感數據的融合方法主要包括簡單或加權平均[10]、主觀分析算法[10]、神經網絡算法[11]和生物光學模型方法[12]等。

經驗正交分解插值方法(Data Interpolating Empirical Orthogonal Function,DINEOF)[13-14]是近年發展起來的一種高效的數據重構方法,它自適應能力強,無需任何先驗知識便可以快速有效地使用的內插方法,在處理高缺失量的大數據集上效率更高,結果更好,在長時序海洋水色遙感數據重構中發揮越來越重要的作用[14-15]。許多國內外學者通過對該方法的不斷改進,滿足不同的重構需求,如Beckers等利用DINEOF和OI方法結合,在完成數據重構的同時,對數據的誤差進行了定量的估計[16];Alvera-Azcárate等利用DINEOF方法結合數據本身的空間相關性對插值數據的“異常值”進行有效剔除[17];Alvera-Azcárate等基于DINEOF方法利用不同變量之間的時空相關性對多變量進行了有效的重構[14]。以上的研究表明,DINEOF方法在海洋水色遙感數據重構中有很強的適用性和改進性。

SeaWiFS和MODIS/Aqua是目前海洋水色遙感研究中廣泛應用的兩種傳感器。SeaWiFS傳感器自1997年9月開始獲取數據,到2010年12月結束壽命,提供了13年的觀測資料;MODIS/Aqua傳感器自2002年7月至今仍在軌運行。兩種數據在2002年~2010年間的同步觀測,有利于兩種數據的對比融合,而兩種數據不同的運行壽命,又為海洋水色參數的長期、連續觀測提供了依據,因此對兩種數據產品的有效重構和對比融合能夠為海洋生態系統的持續研究提供數據支持。

黃、渤海區域(圖1)是中國最重要的半封閉陸架淺海區,其葉綠素a濃度的遙感反演受多種因素的影響,精度和覆蓋率均較低,因此影響了遙感數據在該海區的應用。目前該海區的海洋葉綠素a的遙感研究,主要集中于遙感反演精度的驗證[15]或者直接利用有缺失數據的遙感產品來粗略的分析葉綠素a濃度的動態機制[16-17],對不同傳感器數據缺乏系統的比較,更是較少地涉及到數據完整性重構的系統研究,因而也限制了該海區水色遙感產品的應用。

本文以傳統的DINEOF方法為基礎,以SeaWiFS和MODIS/Aqua海表葉綠素a濃度數據的一致性為依據,借助多變量組合DINEOF方法,對兩種水色遙感數據進行了有效的組合重構,獲得了自1997年9月至2012年6月期間較為完整的海表葉綠素a濃度數據集,為海洋水色遙感數據的應用提供一定的借鑒價值和方法依據。

圖1 研究區

2 數據與方法

2.1 數據來源

本文應用的數據來自于美國國家航空航天局提供的SeaWiFS和MODIS/Aqua(MODISA)海表葉綠素a(Chl-a)遙感產品,考慮到數據覆蓋率和數據處理速度,本文采用全球標準算法8天合成的海表葉綠素a三級數據產品作為實驗數據源。提取兩種數據在重疊時段(2002年7月4日至2010年12月18日) 的所有數據,共388個時段,兩種數據的空間分辨率均為9km×9km,原始數據通過掩膜處理,提取117oE~127oE和31oN~41oN范圍內的數據。

2.2 DINEOF重構的基本原理

DINEOF方法是一種基于經驗正交分解的插值方法,它不需要先驗知識,在處理大數據量和高缺失率數據時有著明顯的優勢[15,18],近年來被越來越多的國外學者用于重構海表溫度[18]、懸浮物質[19]以及葉綠素a[20]等遙感數據。其基本原理是借助經驗正交分解方法,通過對數據集的多次迭代分解和合成,獲得最小的交叉驗證誤差,進而獲得對缺失數據點的最佳重構。具體方法在相關文獻中已經有了詳細的介紹[15-16,18],這里不再詳述。這種方法最主要的優點:①該方法是一種自組織、自適應的數據重構方法,在重構前不需要對數據先驗知識獲取;②該方法基于數據集中所有數據點、所有時段的信息,因此插值結果更能反映數據集的整體特征;③該方法操作簡單,可控性強,與傳統方法相比,運算時間大大縮短(如運算時間僅是最佳內插方法的1/30)。

本文對DINEOF方法重構過程中的相關參數做以下設置:從原始數據集隨機抽取3%的數據點作為交叉驗證數據;EOF分解迭代終止的準則為:前次分解的均方根誤差與本次分解的均方根誤差的差值大于前次均方根誤差的1e-4倍。整個DINEOF重構方法的源代碼基于matlab平臺自行編寫。

2.3 DINEOF組合重構的原理

DINEOF組合重構方法是利用兩組或者多組變量之間的時空相關性,對組合后的數據進行DINEOF重建,提高數據的重構效果,其最初是被用來對海表溫度、葉綠素a濃度和海面風場數據等短期不同變量的數據集進行重構[14]。本文首次嘗試將該方法應用于長時間序列的多源水色遙感相同變量(葉綠素a)產品的組合重構,以彌補單一時段某一數據集的完全缺失現象。組合數據集為兩種數據空間維上的組合:

(1)

其中,XA為組合后的數據集矩陣,XS為SeaWiFS/Chl-a數據集的初始矩陣,XM為MODISA/Chl-a數據集的初始矩陣。

2.4 重構精度的評價

由于大尺度空間的水色遙感數據缺乏有效的實測數據來對重構結果進行評價,當前的精度評價方法多是基于數據集自身進行[19,21-22]。評價方法建立在初始圖像中有數據的象元(有效數據點)基礎上,每次重構之后,計算有效數據點的實際值和重構值之間的相關參數來評價重構效果,這些評價參數包括相關系數、均方根誤差、相對偏差、方差貢獻率等。

圖2 SeaWiFS和MODISA數據的時空覆蓋率

3 結果與討論

黃、渤海海區(圖1)是典型的近岸二類水體,葉綠素a遙感產品受其他懸浮物質的影響,絕對精度較低,需要結合實測資料進行驗證或者建立適合該區域的反演算法[23]。本文承認這種誤差的存在,但是不將其作為研究的內容,而僅側重于對SeaWiFS/Chl-a和MODISA/Chl-a全球算法產品的對比,及DINEOF傳統和組合重構方法對兩種數據重構效果的分析。

3.1 兩種數據時空覆蓋率的比較

圖 2 顯示了兩種數據在研究時段內的時間覆蓋率和空間覆蓋率。圖2(a)、圖2(b)是兩種數據的時間覆蓋率的空間分布特征(圖中的黑色實線是覆蓋率為5%的等值線),兩種數據時間覆蓋率的空間分布形勢基本一致,近岸地區覆蓋率明顯低于離岸地區,最低覆蓋率出現在渤海西南部和江蘇近岸淺灘等水深較淺的區域。圖2(c)、圖2(d)是兩種數據的空間覆蓋率隨時間變化的統計圖,兩種數據空間覆蓋率的時間波動明顯,SeaWiFS數據的空間覆蓋率變化范圍為0.43±0.26(圖2(d)),自2007年之后,空間覆蓋率出現明顯減小趨勢,甚至出現無有效數據覆蓋的時期(2008年~2009年);而MODISA數據(圖2(c))的平均空間覆蓋率為0.49±0.24,略高于SeaWiFS數據,但是波動較小??傮w來看,兩種數據的時空覆蓋率存在著顯著的同步性(R=0.79,p<0.001)。因此,針對這兩種數據產品的融合并不能較大程度的提高數據的時空覆蓋率,但是可以有效地填補單一傳感器(SeaWiFS)某些時段的數據完全缺失狀況,提高其在長時間序列分析中的可用性。

表1 兩種Chl-a數據的數值統計特征

3.2 兩種數據數值特征的比較

將SeaWiFS和MODISA兩種數據的時空重疊部分的數值特征進行統計比較,對兩種數據的一致性進行評價,兩種數據集共有9,311,661個有效重疊象元。表1是兩種數據原始和以10為底的對數轉化后的數值統計特征。結果顯示兩種原始數據集的偏度遠大于0,呈現明顯的正偏態;SeaWiFS數據的峰度偏大,說明數據分布更為集中;MODISA數據的平均值、中值、標準差和變異系數均略高于SeaWiFS數據集,說明MODISA數據集所反演的葉綠素a濃度值略高于SeaWiFS數據集,并且相對變動也較大。進行以10為底的對數轉化后,數據更接近正態分布,呈現微弱的負偏態。因此,兩種數據均呈現對數正態分布[24],在重構和統計分析前,需要對數據進行轉化處理。

圖3 SeaWiFS和MODISA數據的對比,灰度指示了數據點的分布密度

圖3為兩種葉綠素a數據的對數散點圖。相關分析結果顯示,兩種數據線性回歸的斜率為1.04,截距為0.02,相關系數為0.88,兩種數據在研究區域內呈現較好的一致性,兩者的均方根偏差(RMSE=0.17)小于遙感數據本身的平均算法誤差(RMSE=0.22)[25]。圖3顯示該海區葉綠素a濃度多集中于1mg·m-3~10mg·m-3的范圍內,在低值區MODISA數據略低于SeaWiFS數據,在高值區MODISA數據略高于SeaWiFS數據。總體上看,兩種數據在該海區有著高度的一致性,尤其考慮到黃渤海區域大部分海區為二類水體,葉綠素a遙感反演算法本身受多種因素影響,數據反演誤差較大[26-27],但是兩者的數值仍然保持高度一致性,所以SeaWiFS/Chl-a和MODISA/Chl-a數據能相互補充,為該海區葉綠素a濃度的長周期研究提供依據。

3.3 DINEOF重構結果

通過前面的結果可知,兩種數據的時空覆蓋率(圖2)存在明顯的缺失情況,SeaWiFS數據的平均缺失率為59.0%,MODISA數據的平均缺失率為52.9%,總體來看,兩種數據的缺失率較高,并且缺失率隨時間變化較大,甚至出現某些時段圖像的完全缺失情況,導致單一數據集在時空尺度上的不完整性,限制了對長周期規律的探索,因此需要對缺失數據進行有效的重構,以保持數據集在時空尺度上的完整性。首先,利用DINEOF方法對兩種數據進行了單獨重構實驗,然后利用SeaWiFS和MODISA兩種數據的數值特征的一致性,采用DINEOF組合重構方法對兩種數據進行重構,有效地彌補單一數據源的數據完全缺失情況。

在進行重構前,對原始數據進行預處理,對數據值取10為底的對數轉換;然后去除過低時間覆蓋率的數據點,僅保留兩種數據時間覆蓋率均大于5%的數據點,將兩種數據中空間覆蓋率小于5% 的圖像去除;最終所有的數據組成M×N維矩陣(M是空間維,N是時間維),經過處理后SeaWiFS數據組成7347×345的矩陣,MODISA數據組成7347×376的矩陣,組合數據組成14694×388的矩陣,然后利用DINEOF方法分別對3種數據集的缺失值進行重構。

圖4 DINEOF重構過程中交叉驗證點的均方根誤差變化

圖4是DINEOF重構的交叉驗證過程(交叉驗證點),箭頭指示了收斂點及模態數和最小均方根誤差。表2是有效數據點重構結果。從結果來看,兩種數據集的單獨重構和兩者的組合重構均取得較好的重構效果,重構誤差均處在合理的范圍內,均方根誤差小于葉綠素a算法本身的誤差(RMSE=0.22)[25]。SeaWiFS數據集利用35個模態可獲得均方根誤差的收斂,而MODISA數據要利用50個模態才能獲得均方根誤差的收斂,原因可能是MODISA數據相對于SeaWiFS數據的分布較為分散,變異較大;組合數據利用39個模態獲得最佳重構效果。整體上來看,3個重構過程的效果基本一致,均方根誤差(RMSE)為0.1左右,明顯小于數據本身算法誤差(0.22)以及兩種數據之間的偏差(0.17),所以基于DINEOF方法的單獨重構和組合重構均取得了理想的效果。

表2 DINEOF重構結果及驗證

圖5 有效數據點的單獨重構值和組合重構值的對比

圖6 SeaWiFS數據和MODISA數據單獨重構結果

圖7 SeaWiFS數據和MODISA數據組合重構結果

為了定量地對比組合重構和單獨重構的結果,在重構之后,將兩種方法在有效數據點的重構值進行了比較,圖5是SeaWiFS和MODISA數據中有效數據點兩種重構值的對比圖。結果顯示,兩種方法的重構值沒有明顯的偏差,組合重構和單獨重構精度基本一致。

圖6是2003年1月25日至2月1日平均葉綠素a濃度數據的單獨重構結果。從圖上可以看出,在存在較大缺失值的情況下,該方法仍能獲得較好的重構結果,在近岸地區,尤其是長江口附近,存在大量缺失數據的情況下,重構數據仍然較好地反映近岸和河口的高值區,符合前人的研究結果[28-31]。圖7是2009年4月23日至30日平均葉綠素a濃度的組合重構結果。從結果看,雖然SeaWiFS數據在該時間段內完全缺失,但是DINEOF的組合重構方法仍然能很好地對其進行重構,重構后的葉綠素a濃度(圖7(b))與MODISA數據時空特征相似,并且很好地突出了近岸和河口地區的高值,以及春季在黃海中部出現的浮游植物藻華現象[32-33]。

4 結束語

SeaWiFS/Chl-a和MODISA/Chl-a數據產品在黃、渤海區域有著相似的統計特征,兩種數據的時空覆蓋率和數值特征均具有顯著的一致性,進一步驗證了MODIS/Aqua數據可以作為SeaWiFS數據在葉綠素a長時間尺度研究上的有效承接。

利用DINEOF傳統方法和組合重構方法對SeaWiFS/Chl-a數據集、MODISA/Chl-a數據集進行了時空缺失值的重構,雖然兩種數據產品在二類水體中的算法存在較大的誤差,但是DINEOF方法對葉綠素a濃度的重構仍然取得了穩定、理想的效果。另外,對兩種數據的組合數據集的DINEOF重構可以有效地填補單一數據集的時間尺度不連續現象,而且并未給重構結果帶來明顯的額外誤差,重構后的完整數據更有利于該海區生態系統的長期、連續研究需求。

DINEOF方法是近十年發展起來的高效的數據插值重構方法,其在海洋水色數據的插值重構中的應用仍在探索和發展,許多國內外學者也針對其具體應用進行改進。另外,該方法目前主要在海洋環境數據重構中應用,但是根據其原理和效果推斷,該方法亦適用于對陸地、大氣等領域時序數據的重構分析,后續的研究也將嘗試該方面的研究。

參考文獻:

[1] GREGG W W,CONKRIGHT M E.Decadal changes in global ocean chlorophyll[J].Geophysical Research Letters,2002,29(15):1730-1733.

[2] TANG D,KAWAMURA H,LEE M A,et al.Seasonal and spatial distribution of chlorophyll-a concentrations and water conditions in the Gulf of Tonkin,South China Sea[J].Remote Sensing of Environment,2003,85(4):475-483.

[3] GREGG W W,CASEY N W,MCCLAIN C R.Recent trends in global ocean chlorophyll[J].Geophysical Research Letters,2005,32(3):L03606.

[4] VANTREPOTTE V M,LIN F.Inter-annual variations in the SeaWiFS global chlorophyll a concentration (1997-2007)[J].Deep Sea Research Part I:Oceanographic Research Papers,2011,58(4):429-441.

[5] BARNES R A,CLARK D K,ESAIAS W E,et al.Development of a consistent multi-sensor global ocean colour time series[J].International Journal of Remote Sensing,2003,24(20):4047- 4064.

[6] TANG S,DONG Q,LIU F.Climate-driven chlorophyll-a concentration interannual variability in the South China Sea[J].Theoretical and Applied Climatology,2011,103(1-2):229-237.

[7] IIDA T,SAITOH S I.Temporal and spatial variability of chlorophyll concentrations in the Bering Sea using empirical orthogonal function (EOF) analysis of remote sensing data[J].Deep Sea Research Part II:Topical Studies in Oceanography,2007,54(23-26):2657-2671.

[8] REYNOLDS R W,SMITH T M.Improved global sea surface temperature analyses using optimum interpolation[J].Journal of Climate,1994,7(6):929-948.

[9] SCHOELLHAMER D H.Singular spectrum analysis for time series with missing data[J].Geophysical Research Letters,2001,28(16):3187-3190.

[10] POTTIER C,GARCON V,LARNICOL G,et al.Merging SeaWiFS and MODIS/Aqua ocean color data in north and equatorial atlantic using weighted averaging and objective analysis[J].IEEE,2006,44(11):3436-3451.

[11] LOYOLA D,COLDEWEY-EGBERS M.Multi-sensor data merging with stacked neural networks for the creation of satellite long-term climate data records[J].EURASIP Journal on Advances in Signal Processing,2012,2012(1):1-10.

[12] MARITORENA S,DANDON O H F,MANGIN A,et al.Merged satellite ocean color data products using a bio-optical model:Characteristics,benefits and issues[J].Remote Sensing of Environment,2010,114(8):1791-1804.

[13] MILES T N,HE R.Temporal and spatial variability of Chl-a and SST on the South Atlantic Bight:Revisiting with cloud-free reconstructions of MODIS satellite imagery[J].Continental Shelf Research,2010,30(18):1951-1962.

[14] ALVERA-AZCARATE A,BARTH A,BECKERS J M,et al.Multivariate reconstruction of missing data in sea surface temperature,chlorophyll,and wind satellite fields[J].Journal of Geophysical Research,2007,112(C3):C03008.

[15] BECKERS J M,RIXEN M.EOF calculations and data filling from incomplete oceanographic datasets[J].Journal of Atmospheric and Oceanic Technology,2003,20(12):1839-1856.

[16] BECKERS J M,BARTH A,ALVERA A.DINEOF reconstruction of clouded images including error maps.Application to the sea-surface temperature around corsican island[M].European Geosciences Union,2006.

[17] ALVERA-AZC R A,SIRJACOBS D,BARTH A,et al.Outlier detection in satellite data using spatial coherence[J].Remote Sensing of Environment,2012,119:84-91.

[18] ALVERA-AZCARATE A,BARTH A,RIXEN M,et al.Reconstruction of incomplete oceanographic data sets using empirical orthogonal functions:Application to the Adriatic Sea surface temperature[J].Ocean Modelling,2005,9(4):325-346.

[19] NECHAD B,ALVERA-AZCAR TE A,RUDDICK K,et al.Reconstruction of MODIS total suspended matter time series maps by DINEOF and validation with autonomous platform data[J].Ocean Dynamics,2011,61(8):1205-1214.

[20] SHAW P T,CHAO S Y.Surface circulation in the South China Sea[J].Deep Sea Research Part I:Oceanographic Research Papers,1994,41(11-12):1663-1683.

[21] GANZEDO U,ALVERA-AZC RATE A,ESNAOLA G,et al.Reconstruction of sea surface temperature by means of DINEOF:A case study during the fishing season in the Bay of Biscay[J].International Journal of Remote Sensing,2011,32(4):933-950.

[22] SIRJACOBS D,ALVERA-AZC RATE A,BARTH A,et al.Cloud filling of ocean colour and sea surface temperature remote sensing products over the Southern North Sea by the data interpolating empirical orthogonal functions methodology[J].Journal of Sea Research,2011,65(1):114-130.

[23] SISWANTO E,TANG J,YAMAGUCHI H,et al.Empirical ocean-color algorithms to retrieve chlorophyll-a,total suspended matter,and colored dissolved organic matter absorption coefficient in the Yellow and East China Seas[J].Journal of Oceanography,2011,67(5):627-650.

[24] CAMPBELL J W.The lognormal distribution as a model for bio-optical variability in the sea[J].J Geophys Res,1995,100(C7):13237-13254.

[25] O'REILLY J E,MARITORENA S,SIEGEL D,et al.Ocean color chlorophyll a algorithms for SeaWiFS,OC2,and OC4:Version 4[R].SeaWiFS Postlaunch Technical Report Series,Part 3.Greenbelt,Maryland:NASA Goddard Space Flight Center,2000.

[26] SATHYENDRANATH S.Remote sensing of ocean colour in coastal and other optically complex waters[R].Reports of the International Ocean-Colour Coordination Group,No.3.Dartmouth,Canada:IOCCG,2000.

[27] GREGG W W,CASEY N W.Global and regional evaluation of the SeaWiFS chlorophyll data set[J].Remote Sensing of Environment,2004,93(4):463-479.

[28] YAMAGUCHI H,KIM H C,SON Y B,et al.Seasonal and summer interannual variations of SeaWiFS chlorophyll a in the Yellow Sea and East China Sea[J].Progress in Oceanography,2012,105:22-29.

[29] SHI W,WANG M.Satellite views of the Bohai Sea,Yellow Sea,and East China Sea[J].Progress in Oceanography,2012,104:30-45.

[30] 伍玉梅,徐兆禮,崔雪森,等.1997~2007年東海葉綠素a質量濃度的時空變化分析[J].環境科學研究,2008,21(6):137-142.

[31] 叢丕福,牛錚,蒙繼華,等.1998~2003年衛星反演的中國陸架海葉綠素a濃度變化分析[J].海洋環境科學,2006,25(1):30-33.

[32] ZHENG X,WEI H,LI K,et al.Analysis of chlorophyll concentration during the phytoplankton spring bloom in the Yellow Sea based on the MODIS data[G].International Conference on Life System Modeling and Simulation,LSMS,2010:254-261

[33] XUAN J L,ZHOU F,HUANG D J,et al.Physical processes and their role on the spatial and temporal variability of the spring phytoplankton bloom in the central Yellow Sea[J].Acta Ecological Inica,2011,31(1):61-70.

猜你喜歡
方法
中醫特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數學教學改革的方法
河北畫報(2021年2期)2021-05-25 02:07:46
化學反應多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學習方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 一级香蕉视频在线观看| 亚洲一区二区无码视频| 国产福利一区在线| 精品天海翼一区二区| 精品小视频在线观看| 91日本在线观看亚洲精品| 亚洲国产精品美女| 国产丝袜无码精品| 国产无码性爱一区二区三区| a毛片基地免费大全| 精品人妻无码中字系列| 波多野结衣的av一区二区三区| 国内精自视频品线一二区| 亚洲视频欧美不卡| 综合色天天| 四虎在线观看视频高清无码| 在线观看国产网址你懂的| 免费大黄网站在线观看| 国产精品亚洲专区一区| 91娇喘视频| 一级香蕉人体视频| 亚洲手机在线| 色综合天天视频在线观看| 国产91九色在线播放| 精品国产黑色丝袜高跟鞋| 在线人成精品免费视频| 国产真实乱子伦视频播放| 国产丝袜第一页| 国产在线欧美| 全部免费特黄特色大片视频| 又粗又大又爽又紧免费视频| 欧美午夜小视频| 免费jjzz在在线播放国产| 伊人成人在线视频| 日韩无码一二三区| 国产日韩av在线播放| 国产无码高清视频不卡| 国产欧美日本在线观看| 国产特级毛片aaaaaaa高清| 99精品伊人久久久大香线蕉| 国产精品极品美女自在线看免费一区二区| 久久久91人妻无码精品蜜桃HD | 97在线国产视频| 亚洲一区二区三区国产精品| 日本免费新一区视频| 国产精品亚洲一区二区在线观看| 在线观看精品国产入口| 国产无码网站在线观看| 国产精品香蕉在线| 亚洲欧美人成人让影院| 欧洲高清无码在线| 99人妻碰碰碰久久久久禁片| 欧美在线视频不卡第一页| 在线播放真实国产乱子伦| 国产精品乱偷免费视频| 日韩av电影一区二区三区四区 | 视频一区亚洲| 亚洲伦理一区二区| 亚洲成人精品| 婷婷亚洲视频| 久久伊伊香蕉综合精品| 色播五月婷婷| 日韩成人高清无码| 中文字幕色站| 网久久综合| 天天综合网亚洲网站| 国产成人高清精品免费| 国产丰满成熟女性性满足视频| 亚洲一级色| 国产精品专区第1页| 香蕉网久久| 真实国产乱子伦高清| 激情综合激情| 中文字幕乱码中文乱码51精品| 欧美日韩国产系列在线观看| 91久久青青草原精品国产| 99热6这里只有精品| 手机在线免费不卡一区二| 午夜福利网址| 日本久久网站| 伊人中文网| 精品一區二區久久久久久久網站|