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

利用環(huán)境一號衛(wèi)星熱紅外影像反演渤海海表溫度

2011-01-09 05:22:56余曉磊巫兆聰
海洋技術(shù)學(xué)報 2011年2期
關(guān)鍵詞:大氣環(huán)境

余曉磊,巫兆聰

(武漢大學(xué) 遙感信息工程學(xué)院,湖北 武漢 430079)

利用環(huán)境一號衛(wèi)星熱紅外影像反演渤海海表溫度

余曉磊,巫兆聰

(武漢大學(xué) 遙感信息工程學(xué)院,湖北 武漢 430079)

針對環(huán)境衛(wèi)星熱紅外遙感影像,結(jié)合美國環(huán)境預(yù)報中心(NCEP)再分析數(shù)據(jù),運(yùn)用覃志豪單窗算法(MW),修訂了該算法的主要參數(shù)計算公式,建立了反演海洋表面溫度的流程。利用2009年10月4日渤海上空的熱紅外遙感影像進(jìn)行反演試驗(yàn),同時對比反演結(jié)果和美國中等分辨率成像光譜儀(MODIS)的海表溫度產(chǎn)品(MOD28)之間的差異。結(jié)果表明,反演得到海表溫度的空間分布趨近一致,與MOD28相比相對誤差在5%左右,具有良好的相關(guān)性。

環(huán)境一號衛(wèi)星;海表溫度反演;單通道算法;MODTRAN;MODIS;MOD28

海洋表面溫度(SST)是研究海洋表面海-空水汽和熱量交換的一個重要物理參量,它同時還影響著海洋生物的生長、繁殖和分布特征以及各類的海洋環(huán)境要素[1]。渤海作為我國的內(nèi)海,具有豐富優(yōu)質(zhì)的漁業(yè)、港口、石油、景觀和海鹽資源,其水溫變化受北方大陸性氣候影響。隨著海洋資源的開發(fā)利用活動,渤海的資源和生態(tài)環(huán)境受到了較大的破壞。渤海環(huán)境質(zhì)量不斷惡化,主要表現(xiàn)于海岸帶污染明顯、污染范圍擴(kuò)大、生態(tài)系統(tǒng)弱化、生態(tài)環(huán)境退化、赤潮、富營養(yǎng)化等[2]。因此,監(jiān)測渤海海面溫度成為迫切必要的任務(wù)。傳統(tǒng)實(shí)地測量的方法無法滿足大面積實(shí)時監(jiān)測海溫的需要,遙感技術(shù)可以大面積重復(fù)觀測海洋表面,為反演海表溫度提供了便利[3,13]。

環(huán)境一號衛(wèi)星(HJ-1)是我國2008年9月6日發(fā)射的專門用于環(huán)境和災(zāi)害監(jiān)測預(yù)報的小衛(wèi)星星座,它由兩顆衛(wèi)星組成(HJ-1A和HJ-1B)。其中HJ-1B搭載的紅外相機(jī)(IRS)可以完成對地幅寬為720 km、地面像元分辨率為150 m/300 m、近短中長4個光譜譜段的成像,重復(fù)周期4 d(表1)。其熱紅外波段與Landsat-TM6號波段相仿。

本文利用IRS影像的熱紅外波段(IRS8),采用覃志豪單窗算法(MW)反演渤海海表溫度,修正了主要參數(shù),針對反演所需的大氣平均作用溫度以及大氣柱狀水汽含量不易實(shí)時獲取的問題,嘗試采用美國環(huán)境預(yù)報中心(NCEP)再分析數(shù)據(jù)進(jìn)行對偶克里金(Kriging)插值。同時將反演結(jié)果和中等分辨率成像光譜儀(MODIS)的 MOD28(SST)產(chǎn)品,進(jìn)行了比較。結(jié)果表明,反演的海表溫度同MOD28相比,相對誤差5%左右,而利用環(huán)境星影像反演溫度的動態(tài)范圍比MODIS影像反演的要大,且其空間分布趨近一致。說明利用環(huán)境衛(wèi)星監(jiān)測海表溫度具有廣闊的應(yīng)用前景。

表1 HJ-1B有效載荷參數(shù)

1 反演算法

1.1 單窗算法及其修訂

對于遙感影像熱紅外波段(10.5~12.5 μm),假設(shè)海洋表面和大氣對熱輻射的傳導(dǎo)有朗伯體特性,則按照輻射傳輸方程,傳感器接收到的熱紅外輻射亮度可以表示為[4]:

式中:Bsensor,λ表示傳感器接收到的熱紅外輻射亮度;Batm↑表示大氣上行輻射亮度;Latm↓表示大氣下行輻射亮度;ελ表示海表在熱紅外波段的發(fā)射率,一般在無浪的情況下可認(rèn)為是 0.98[5,13];τλ表示大氣在熱紅外波段的大氣透射率;Bλ(Ts)表示海表自身的輻射亮度,其中Ts表示海表溫度。同時根據(jù)輻亮度定義以及普朗克公式有[5]:

式中:c為光速;h為普朗克常數(shù);k為波爾茲曼常數(shù),c1=2πhc2=3.741 8×10-16W·m2,c2=hc/k=1.438 8×104μm·K。

考慮大氣上行輻射亮度與下行輻射亮度大致相等的條件下,式(1)可以簡化為[6]:

式中:Ta表示大氣平均作用溫度。

覃志豪等通過對輻射傳輸方程線性展開,針對TM6影像,建立了適用于熱紅外通道(10.5~12.5 μm)反演溫度的單窗算法(Mono-Window Algorithm)[6-7]:

式中:C=ελ·τλ,D=(1-τλ)[1+τλ(1-ελ)],覃志豪等在 MW算法中對普朗克公式(式3)在地球表面的溫度范圍內(nèi)(0℃~70℃)按照泰勒級數(shù)線性展開,并據(jù)此定義溫度參數(shù)L:

發(fā)現(xiàn)參數(shù)L與溫度T有很好的線性關(guān)系,根據(jù)這一關(guān)系建立L于T之間的聯(lián)系(如式7),對于HJ-1B衛(wèi)星的IRS8波段的 a,b 為回歸系數(shù),依據(jù)式(3),(6),(7),當(dāng)海表溫度在-2℃~70℃之間時,a=-68.690 7,b=0.586 6,如圖 1 所示,相關(guān)系數(shù)的平方為0.999 8,RMSE=0.223 4。

圖1 參數(shù)L隨溫度T變化關(guān)系圖

對于大氣平均作用溫度Ta,利用MODTRAN4輻射傳輸模型中自帶的美國1976標(biāo)準(zhǔn)大氣、熱帶大氣、中緯度夏季大氣和中緯度冬季大氣的大氣剖面參數(shù),并根據(jù)Ta的近似計算公式[6-7,15]:

式中:m為大氣剖面層數(shù);Tz為每層大氣溫度;w為水汽總含量;w(z)為每層大氣水汽含量;Rw(z)為每層水汽占總含量之比,且

其中:T0為近地表溫度;Rt(z)為逐層大氣溫度隨高度變化的比率。依照式(8)~式(9)并集合常見大氣剖面數(shù)據(jù)可得表2,以計算大氣平均作用溫度Ta。

表2 IRS8通道大氣平均作用溫度計算表

對于熱紅外波段的大氣透射率τλ,通過分析整理大連探空站 (38°57′43″N~121°25′05″N)2008 年 9 月 1 日—2008 年10月31日的大氣探空數(shù)據(jù),運(yùn)用大氣輻射傳輸計算程序MODTRAN4結(jié)合IRS8通道的光譜響應(yīng)函數(shù)(如圖2)模擬其同大氣柱狀水汽含量w的關(guān)系后,有以下線性近似關(guān)系[7,15]:

相關(guān)系數(shù)的平方為0.976 8,RMSE=0.003 4。

1.2 NCEP再分析數(shù)據(jù)的Kriging插值

如前所述,用單窗算法反演海表溫度,需要獲取近地表的大氣溫度和大氣柱狀水蒸氣含量數(shù)據(jù)。對于海洋地區(qū),其上空的無線電探空資料有限,無法滿足實(shí)時反演的需要。美國國家環(huán)境預(yù)報中心(NCEP)和國家大氣研究中心(NCAR)聯(lián)合執(zhí)行的全球大氣資料再分析計劃,通過收集全球無線電測風(fēng)資料、綜合海氣資料(COADS)、飛機(jī)觀測資料、陸面天氣觀測資料、衛(wèi)星探測資料、微波特殊探測/圖像(SSM/I)資料和衛(wèi)星觀測風(fēng)(云的移動)資料等一系列數(shù)據(jù),利用巨型計算機(jī)進(jìn)行數(shù)據(jù)同化(Data Assimilation),并對外發(fā)布每日逐6 h再分析資料集[8]。數(shù)據(jù)集為全球1°×1°經(jīng)緯格網(wǎng)點(diǎn)處的位勢高度、溫度廓線、濕度廓線、氣壓廓線和太陽輻射通量等數(shù)據(jù)。

本文利用NCEP分析數(shù)據(jù),獲取渤海上空4個格網(wǎng)點(diǎn)的近地表溫度和大氣柱狀水汽含量,同時采用對偶Kriging插值[9],以獲取每個象元對應(yīng)近地表溫度和大氣柱狀水蒸汽含量值。

2 算法實(shí)現(xiàn)

2.1 環(huán)境衛(wèi)星數(shù)據(jù)預(yù)處理

通過中國資源衛(wèi)星應(yīng)用中心獲取了2009年10月4日過境的HJ-1B紅外相機(jī)(IRS)二級影像,掃描時間為 5:34(UTC) 截取了覆蓋范圍為 37°06′08″N~40°11′08″N,117°45′04″E~121°53′49″E 的渤海上空的子影像。以 1:50 000 的基礎(chǔ)地形圖為基準(zhǔn)采用三次多項(xiàng)式模型,雙線性內(nèi)插采樣方法對原始衛(wèi)星影像進(jìn)行了幾何糾正,糾正后的圖像以WGS-84橢球?yàn)閰⒖迹队胺绞綖閁TM投影,糾正精度在0.3個像元以內(nèi)。將幾何糾正后的影像進(jìn)行云掩模和海陸淹沒運(yùn)算,只保留海表像元。

在反演海表溫度之前需要把影像中每個像元的數(shù)字灰度值(DN)轉(zhuǎn)換為亮度溫度。首先根據(jù)資源衛(wèi)星應(yīng)用中心提供的定標(biāo)公式和定標(biāo)系數(shù),將DN值轉(zhuǎn)換為輻射能量強(qiáng)度,定標(biāo)公式:

對于HJ-1B的IRS8波段,式中b=-25.441為偏移量,a=59.421(DN/W·m-2·sr-1·μm-1)為增益,DN 為數(shù)字灰度值。

根據(jù)式(3)可以推導(dǎo)得:

式中:K1=2hc2/λ5(W·m-2·sr-1·μm-1),K2=hc/kλ(K),其中λ為IRS8波段的等效中心波長,按資源衛(wèi)星應(yīng)用中心提供的 IRS8 波段響應(yīng)函數(shù)(如圖 2),同時依照式(13)[5],計算等效中心波長:

其中 f(λ)為 IRS6波段的光譜響應(yīng)函數(shù),λmax與 λmin為光譜波段范圍,由此可以計算出λc=11.494 8 μm,則:

2.2 海溫(SST)反演流程

海溫反演流程如圖3所示:

圖2 HJ-1B IRS8波段光譜響應(yīng)曲線

圖3 HJ-1B IRS8波段反演海表溫度(SST)流程

3 結(jié)果分析

3.1 中等分辨率成像光譜儀(MODIS)反演海表溫度算法簡介

中等分辨率成像光譜儀(MODIS)作為美國對地觀測系統(tǒng)(EOS)中的主要傳感器搭載在TERRA和AQUAR衛(wèi)星上,可以每天兩次獲取全球任意地點(diǎn)的影像數(shù)據(jù),其中的31號(中心波長 11.03 μm),32 號波段 (中心波長 12.202 μm),類似于NOVAA衛(wèi)星AVHRR傳感器的4、5號波段,處于大氣的熱紅外窗口,受到太陽光反射的影響微弱,且對于水汽的吸收作用不同,可以用來校正水汽吸收的影響,其反演海表溫度(SST)采用分裂窗口算法[10]。NASA下屬的海洋水色工作組,將海表溫度(SST)作為 MODIS的二級產(chǎn)品(MOD28)對外分發(fā)[11],經(jīng)過長期的大洋浮標(biāo)數(shù)據(jù)以及船舶報數(shù)據(jù)的驗(yàn)證表明,其精度為 0.053℃~0.66℃[10-11]。

3.2 環(huán)境星反演結(jié)果和MOD28的比較

圖4 HJ-1-B反演SST結(jié)果(a)和MODIS反演SST結(jié)果(b)對比

表3 環(huán)境星反演結(jié)果與MOD28比較表

由于缺少海表實(shí)測溫度數(shù)據(jù),所以將用HJ-1-B衛(wèi)星IRS影像反演得到的海表溫度和MODIS的標(biāo)準(zhǔn)海表溫度(如圖 4)產(chǎn)品(MOD28)進(jìn)行交叉驗(yàn)證[14](cross validation),隨機(jī)抽取30個點(diǎn),比較二者之間的差異,如表3所示。

其中MOD28海表溫度產(chǎn)品由2009年10月4日03:22: 02(UTC)過境的 TERRA 衛(wèi)星 MODIS傳感器獲取,與HJ-1B 衛(wèi)星同日 02:52:31(UTC)的過境時間相差約 0.5 h,基本可以認(rèn)為兩時刻海表溫度保持不變。

圖5 環(huán)境星反演結(jié)果與MOD28之間的相關(guān)性

從表3可以看出,環(huán)境衛(wèi)星反演的海表溫度與MOD28相比,最大偏差在5%以內(nèi),由圖4可以發(fā)現(xiàn),環(huán)境衛(wèi)星反演結(jié)果在15.0℃~27.2℃之間,動態(tài)范圍大于MOD28(15.1℃~26.5℃之間)。進(jìn)一步分析二者之間的相關(guān)性,如圖5,其相關(guān)性為0.936 2,標(biāo)準(zhǔn)偏差為1.374。可見環(huán)境衛(wèi)星反演的海表溫度與MOD28有較高的一致性。

3.3 渤海海域海表溫度分布特征

通過圖4可發(fā)現(xiàn),渤海海域在2009年10月4日海表溫度呈南高北低的總體趨勢,而且溫度分布趨勢沿岸低,中心區(qū)域高,除有明顯人為活動作用的沿岸港口外,中心呈閉合暖中心區(qū)域的溫度比沿岸海域的溫度高0.68-2.74℃。但是遼東半島老鐵山附近海域有一明顯的冷舌區(qū),整個地區(qū)最低溫度達(dá)16.34℃。反演結(jié)果的海表溫度分布規(guī)律同賈瑞麗,孫璐[12]對黃、渤海冬夏季主要月份的海溫分布特征相符,表明衛(wèi)星反演海表溫度分布和多年觀測平均值具有一致性。

如圖6所示,分別截取了天津塘沽港(a)以及唐山市曹妃甸新區(qū)(b)的環(huán)境衛(wèi)星反演海表溫度影像。根據(jù)影像統(tǒng)計,塘沽港附近海面的溫度平均比外海區(qū)域海表溫度高3.69℃,這主要是由頻繁往來的船只,船艙冷卻水對外排放以及港口地區(qū)生產(chǎn)生活產(chǎn)生的熱水所致[2]。而唐山市曹妃甸工業(yè)區(qū),由于近年來填海造陸,新建鋼鐵廠和港口,人類活動頻繁,其附近海域的海表溫度平均比外海區(qū)域高4.01℃。可見渤海海域沿岸由于人類活動加巨,工業(yè)生活快速發(fā)展,對海表溫度的影響巨大,從而間接的影響到該海域,特別是近海岸地帶的海洋生態(tài)環(huán)境。

圖6 環(huán)境星反演塘沽口(a)與曹妃甸(b)海域海表溫度

3.4 誤差分析

對于單窗算法反演海表溫度(SST),反演精度受限于大氣透射率和環(huán)境溫度的修正,這需要掌握反演區(qū)域上空實(shí)時的大氣狀態(tài)廓線。實(shí)際對于廣闊的海洋區(qū)域,缺乏探空資料。本文嘗試采用美國環(huán)境預(yù)報中心(NCEP)再分析數(shù)據(jù),通過對偶克里金插值,獲取反演區(qū)域的大氣狀態(tài)資料,以確定反演所需要的大氣平均作用溫度和水汽含量參數(shù),但是由于所覆蓋的NCEP再分析數(shù)據(jù)只有1°×1°分辨率,對大范圍的影像仍顯稀少,所以直接影響了反演的精度。在本文討論的反演流程中,認(rèn)為海表的發(fā)射率ε為0.98,但是海水的發(fā)射率還可能隨水泥沙含量、海浪狀況,和觀測幾何條件[5]而發(fā)生改變。相應(yīng)的研究表明,高視角天頂角條件下,平靜海面的發(fā)射率可能下降到0.95,這也影響了反演的精度。此外,由于海洋表面的“皮膚效應(yīng)”,衛(wèi)星遙感反演的海表溫度實(shí)際上是海面幾微米薄層內(nèi)的平均溫度,而傳統(tǒng)的海溫測量,則是量測海面0.1~1 m深的水層的溫度[3,5],在實(shí)際運(yùn)用過程中,還需要作出相應(yīng)的修正。

4 結(jié)論與展望

(1)利用環(huán)境衛(wèi)星熱紅外遙感影像,修訂了基于TM6波段的單窗算法,反演海表溫度,反演結(jié)果具有良好的時空一致性,同時與MOD28相比,具有很高的相關(guān)性。

(2)使用NCEP再分析數(shù)據(jù)進(jìn)行克里金插值獲取研究區(qū)域上空的大氣狀況可用于確定單窗算法反演海表溫度所需的參數(shù),對于進(jìn)一步業(yè)務(wù)化流程的實(shí)現(xiàn)具有推廣意義。

(3)由于缺乏實(shí)測海表溫度資料,本文僅對環(huán)境衛(wèi)星反演結(jié)果和MOD28產(chǎn)品進(jìn)行了交叉驗(yàn)證,在下一步的研究應(yīng)用中,還需結(jié)合實(shí)測資料,對算法的精度和可靠性作出進(jìn)一步的分析和評價。

[1]SEELYE M.Introduction to ocean remote sensing[M].Cambridge,UK:Cambridge University Press,2004.

[2]劉學(xué)海.渤海近岸水域環(huán)境污染狀況分析[J].環(huán)境保護(hù)科學(xué),2010,36(1):14-18.

[3]劉良明.衛(wèi)星海洋遙感導(dǎo)論[M].武漢:武漢大學(xué)出版社,2005.

[4]毛克彪,唐華俊,周清波,等.用輻射傳輸方程從MODIS數(shù)據(jù)中反演地表溫度的方法[J].蘭州大學(xué)學(xué)報(自然科學(xué)版),2007,43(4):12-17.

[5]田國良,等.熱紅外遙感[M].北京:電子工業(yè)出版社,2006.

[6]覃志豪,張明華,ARNON K,等.用陸地衛(wèi)星TM6數(shù)據(jù)演算地表溫度的單窗算法[J].地理學(xué)報,2001,56(4):456-466.

[7]覃志豪,李文娟,張明華,等.單窗算法的大氣參數(shù)估計方法[J].國土資源遙感,2003,2:37-43.

[8]KALNAY E,KANAMITSU M,KISTLER R,et al.The NCEP/NCAR 40-year reanalysis project[J].Bulletion of the American Meteorological Society,1996,77,437-470.

[9]鄭永駿,金之雁.對偶Kriging插值方法在氣象資料分析中的作用[J].應(yīng)用氣象學(xué)報,2008,9(2):201-208.

[10]MODIS遙感信息處理原理與算法[M].北京:科學(xué)出版社,2001.

[11]OTIS B B,PETER J M,et al.Modis Infrared Sea Surface Temperature Algorithm Theoretical Basis Document[EB/OL].[2004-5-6].http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=28

[12]賈瑞麗,孫璐.渤海、黃海冬夏季主要月份的海溫分布特征[J],海洋通報,2002,21(4):1-8.

[13]于杰,李永振,陳丕茂,等.利用LANDSAT TM6數(shù)據(jù)反演大亞灣海水表層溫度[J].國土資源遙感,2009,3:24-29.

[14]劉恒.多傳感器衛(wèi)星海表溫度數(shù)據(jù)的驗(yàn)證與交叉比較[D].青島:中國海洋大學(xué),2008.

[15]毛克彪,覃志豪.大氣輻射傳輸模型及MODTRAN中透射率計算[J],測繪與空間地理信息,2004,27(4):1-2.

SST Retrieving of Bohai Sea Using Thermal Infrared Images on HJ-1 Satellite

YU Xiao-lei,WU Zhao-cong
(School of Remote Sensing and Information Engineering,Wuhan University,Wuhan Hubei 430079,China)

Using thermal infrared images from HJ-1 Satellite,combined with the National Centers for Environmental Prediction(NCEP)reanalysis data,the sea surface temperature(SST)retrieving process was established by mono-window algorithm(MW)with modifying the main parameter calculating formulas.The thermal infrared remote sensing image obtained on October 4th 2009,over the Bohai Sea was utilized to validate the process.The comparison between the result made by HJ-1 Satellite and the Sea Surface Temperature product(MOD28)made by MODIS showed that the distribution characteristics of SST is consistent to the former research and the most warps between them is less than 5%,which indicates that they have a nice correlation.

HJ-1 Satellite;retrieval of SST;mono-window algorithm;MODTRAN;MODIS;MOD28

TP79

A

1003-2029(2011)02-0001-06

2011-04-15

國家高技術(shù)研究發(fā)展計劃(863計劃)資助項(xiàng)目(2007AA120203):遙感軟件體系架構(gòu)及標(biāo)準(zhǔn)規(guī)范研究。

余曉磊(1986-),男,湖北武漢人,博士研究生,主要從事遙感信息定量化方面研究。

致謝:感謝中國資源衛(wèi)星應(yīng)用中心提供的環(huán)境衛(wèi)星影像數(shù)據(jù)、定標(biāo)數(shù)據(jù)、光譜響應(yīng)函數(shù)文件。

猜你喜歡
大氣環(huán)境
大氣的呵護(hù)
軍事文摘(2023年10期)2023-06-09 09:15:06
太赫茲大氣臨邊探測儀遙感中高層大氣風(fēng)仿真
長期鍛煉創(chuàng)造體內(nèi)抑癌環(huán)境
一種用于自主學(xué)習(xí)的虛擬仿真環(huán)境
孕期遠(yuǎn)離容易致畸的環(huán)境
不能改變環(huán)境,那就改變心境
環(huán)境
孕期遠(yuǎn)離容易致畸的環(huán)境
大氣古樸揮灑自如
大氣、水之后,土十條來了
主站蜘蛛池模板: 国产欧美日韩91| 亚洲欧美精品一中文字幕| 国产欧美在线观看一区| 丰满人妻中出白浆| 国产爽歪歪免费视频在线观看| 91精品国产无线乱码在线 | 国产精品露脸视频| 国产性精品| 91网站国产| 婷婷六月综合| 幺女国产一级毛片| 乱人伦视频中文字幕在线| 精品国产女同疯狂摩擦2| 精品欧美视频| 国产精品蜜臀| 亚洲男人的天堂在线观看| 久久综合九色综合97婷婷| 一级毛片在线直接观看| 久久精品人人做人人爽97| 国产91麻豆免费观看| 色香蕉网站| 久久国产免费观看| 色哟哟精品无码网站在线播放视频| 97国内精品久久久久不卡| 福利小视频在线播放| 丝袜亚洲综合| 欧美在线一二区| 午夜不卡视频| 亚洲成a人片在线观看88| 中国一级毛片免费观看| 成人免费视频一区| 亚洲国产综合精品中文第一| 亚洲成a人片77777在线播放| 成人噜噜噜视频在线观看| 女人毛片a级大学毛片免费| 中文天堂在线视频| 丁香婷婷激情网| 国产成人免费| 欧美日韩中文国产| 成人午夜视频在线| 一级成人a做片免费| 99re在线观看视频| 成人一级免费视频| 日韩第一页在线| 2020国产精品视频| 久久精品国产91久久综合麻豆自制| 波多野结衣一区二区三区AV| 国产精品欧美激情| 欧美性色综合网| 国产精品粉嫩| 丁香五月婷婷激情基地| 精品伊人久久久久7777人| 亚洲午夜国产片在线观看| 国产香蕉国产精品偷在线观看 | 国产在线视频二区| 最新国产午夜精品视频成人| 欧美一区二区三区不卡免费| 日韩欧美视频第一区在线观看| 98精品全国免费观看视频| 日韩无码白| 日韩小视频网站hq| 有专无码视频| 幺女国产一级毛片| 色噜噜狠狠狠综合曰曰曰| 国产精品性| 91亚洲视频下载| 自拍偷拍欧美日韩| 欧美色图久久| 亚洲无码高清视频在线观看 | 夜夜操国产| 午夜国产大片免费观看| 扒开粉嫩的小缝隙喷白浆视频| 午夜在线不卡| 免费aa毛片| 亚洲成a人片在线观看88| 怡春院欧美一区二区三区免费| 中文无码影院| 免费无码AV片在线观看中文| 热re99久久精品国99热| 日韩欧美一区在线观看| 亚洲av成人无码网站在线观看| 久久婷婷六月|