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

中紅外大氣強吸收通道的地表輻亮度圖像模擬

2017-09-21 01:11:16劉瑤張文娟張兵甘甫平
自然資源遙感 2017年3期
關鍵詞:大氣模型

劉瑤, 張文娟, 張兵, 甘甫平

(1.中國國土資源航空物探遙感中心,北京 100083; 2.中國科學院遙感與數字地球研究所數字地球重點實驗室,北京 100094)

中紅外大氣強吸收通道的地表輻亮度圖像模擬

劉瑤1, 張文娟2, 張兵2, 甘甫平1

(1.中國國土資源航空物探遙感中心,北京 100083; 2.中國科學院遙感與數字地球研究所數字地球重點實驗室,北京 100094)

針對SPIRIT-Ⅲ遙感器的2個4.3 μm強吸收通道,以鄰近4.3 μm的MODIS第23通道圖像為數據源,進行地表輻亮度圖像模擬。首先,通過建立通道轉換模型生成強吸收通道的地表發射率圖像; 然后,基于已有的中紅外輻射傳輸解析模型,推導獲得4.3 μm吸收通道的地表輻亮度解析模型,基于地表的解析模型和輻射傳輸模型MODTRAN的模擬數據,提出吸收通道的大氣效應參數計算方法; 最后,基于地表輻亮度解析模型,綜合發射率和溫度圖像、大氣效應參數、以及通道響應函數實現2個吸收通道的地表輻亮度圖像模擬。對模擬數據進行精度評價,結果表明地表發射率模擬誤差不超過±6%,地表輻亮度模擬誤差不超過±0.02%。由此可見,本文提出的4.3 μm強吸收通道的地表輻亮度圖像模擬方法可行,能夠為該吸收通道的入瞳輻亮度圖像和遙感器成像過程模擬提供地表輻射數據支撐。

圖像模擬; 中紅外; 大氣吸收; 地表發射率

0 引言

以4.3 μm為中心吸收波長的中紅外大氣強吸收波段,主要應用在紅外預警系統中,用于導彈目標的探測與跟蹤[1]。早在20世紀80年代,美國的國防支援計劃(defense support program,DSP)已經開始將該吸收波段作為其紅外傳感器的工作波段。目前,4.3 μm的大氣強吸收通道仍應用在取代DSP的天基紅外系統(space based infrared systems,SBIRS)當中[2]。SBIRS系統已分別于2011年和2013年發射2顆衛星,預計2016年完成全部4顆衛星的發射。由此可見,4.3 μm強吸收通道的傳感器研發和制造對于預警系統的建設十分重要。圖像模擬已成為新型遙感器設計和優化的重要手段[3-6]。4.3 μm強吸收通道的圖像難以獲取,利用圖像模擬手段獲得研究數據是較為可行的方法之一。但是,目前國內外的圖像模擬模型和軟件主要集中在可見光-短波紅外以及熱紅外波段,中紅外波段圖像模擬技術研究相對較少。遙感圖像模擬分為地表場景模擬、大氣模擬和傳感器成像過程模擬3個環節。考慮到傳感器成像過程模擬與傳感器的載荷參數密切相關,在無法獲得傳感器詳細信息的條件下,本文將研究重點放在地表場景和大氣過程的模擬方面。

本文選擇空間紅外成像望遠鏡(spatial infrared imaging telescope Ⅲ,SPIRIT-Ⅲ[7])的2個4.3 μm強吸收通道S1(4.21~4.37 μm)和S2(4.23~4.47 μm)為研究對象,提出一種基于MODIS第23通道(以下簡稱MODIS23通道)(4.02~4.08 μm)數據的輻亮度圖像模擬方法。SPIRIT-Ⅲ是1996年發射的中段空間試驗衛星(midcourse space experiment,MSX)上搭載的紅外遙感器。MSX的主要任務是為導彈防御工作收集目標和背景數據[8]。根據已有研究,4.3 μm的大氣強吸收通道地表輻亮度需同時考慮地表反射與發射,且地表輻亮度與地表光譜特性(發射率或發射率)、地表溫度、大氣狀況及觀測幾何相關[9]。鑒于溫度圖像可通過溫度反演算法獲取[10],觀測幾何數據亦有計算方法[11],本文研究旨在解決4.3 μm強吸收通道的地表發射率和大氣效應的模擬問題,以實現該通道的地表輻亮度圖像模擬,從而為強吸收通道的入瞳輻亮度模擬、遙感器成像過程模擬等后續工作提供數據支撐。

1 數據與方法

1.1 4.3 μm強吸收通道的地表發射率模擬

本文基于地物發射率在鄰近通道間的相關性[12],建立4.3 μm鄰近的非吸收通道數據(即MODIS23通道)到待模擬吸收通道的地表發射率通道轉換模型,從而實現4.3 μm強吸收通道的地表發射率圖像模擬。用于建立發射率通道轉換模型的樣本數據來自加州大學圣塔芭芭拉分校的MODIS發射率光譜庫(MODIS UCSB Emissivity Library)[13]。

選取MODIS UCSB發射率光譜庫中的137種常見類型(植被、水體、土壤和人工類別)光譜,計算它們在MODIS23通道和SPIRIT通道的等效發射率。對于通道i,等效發射率ε(i)的計算方法為

(1)

式中:R(λ)和ε(λ)分別為波長λ處的光譜響應函數值和地物發射率;λi,L和λi,U為通道的起始和終止波長。MODIS和SPIRIT通道[14]的光譜響應函數如圖1所示。

圖1MODIS23通道和SPIRIT-ⅢS1和S2通道的響應函數

Fig.1SpectralresponsefunctionsofMODIS23bandandSPIRIT-ⅢS1,S2bands

根據公式(1)計算得到MODIS23通道和SPIRIT-Ⅲ兩個吸收波段的通道等效發射率,結果如圖2所示。由圖2可見,地物在MODIS23和S1,S2通道的發射率非常接近,說明數據源通道與模擬通道的發射率之間應具有很好的相關性。

圖2 各類地物在MODIS23通道和模擬通道的等效發射率

選擇其中114種地物作為樣本數據,用于分析MODIS23通道發射率和S1/S2通道發射率的相關關系,構建發射率轉換模型,結果如圖3所示。

圖3 MODIS的23通道發射率和S1/S2通道的發射率相關性

由圖3可知,從MODIS23通道向S1/S2通道的等效發射率轉換模型分別為

εS1=0.915 1εMODIS23+0.081 9,

(2)

εS2=0.816 8εMODIS23+0.176 1,

(3)

式中: MODIS23通道的等效發射率為εMODIS23; S1和S2通道的等效發射率分別為εS1和εS2。圖3還給出了通道轉換模型的決定系數R2。2個轉換模型的決定系數都大于0.97,表明MODIS23和S1/S2通道的等效發射率具有很好的線性相關。因此,S1/S2通道的發射率圖像,可基于MODIS23的地表發射率圖像,采用公式(2)(3)求得。

1.2 4.3 μm強吸收通道的地表輻亮度模擬模型

1.2.1 地表輻亮度的輻射傳輸解析模型

Cota等于2009年提出了中紅外輻射傳輸解析模型[15],并將該模型成功應用在參數化成像鏈分析與模擬軟件(Parameterized Image Chain Analysis & Simulation Software,PICASSO)中[16]。該輻射傳輸模型考慮了目標像元和背景像元的差異,且適用于大氣強吸收波段,本文將通過改進這一模型以獲得吸收通道的地表輻亮度解析模型。

根據原始的輻射傳輸解析模型[15],在已知地表信息(溫度、發射率/反射率)和表征大氣效應的相關參數的前提下,能夠模擬計算中紅外波長λ的入瞳輻亮度。對于地表輻亮度,由于不包括大氣上行的輻射傳輸過程,可理解為傳感器位于地表時的入瞳輻亮度。因此,由背景像元散射后直接到達入瞳處的輻亮度和大氣程輻射項都無需考慮。由此可知,地表輻亮度只包含原輻射傳輸模型中入瞳輻亮度L(λ)的第1,4,7,9這4項,即

(4)

(5)

式中:LBB(λ,Tt)是波長λ和地表溫度Tt時的Planck函數,其中h為普朗克常數;c為光速;k為玻爾茲曼常數;ρt(λ)和ρb(λ)分別為目標像元和背景像元在波長λ的反射率;εt(λ)和εb(λ)是對應的發射率;Tt和Tb是對應的地表溫度;A(λ)和A′(λ)分別是目標像元反射的太陽輻射和大氣熱輻射;S(λ)為地表的大氣半球反照率;τd(λ)為地表到大氣頂的上行直射透過率; 其他4個參數是表征大氣效應的物理量,通常稱之為大氣效應參數,指的是只取決于輻射源(一般指太陽)、大氣狀況和觀測幾何,而與地表條件無關的變量。

地表輻亮度不受到上行直射透過率τd引起的衰減。故地表輻亮度LBOA(λ)的表達式為

(6)

(7)

考慮實際計算往往不區分目標和背景,假設ρt(λ)=ρb(λ)=ρ(λ),εt(λ)=εb(λ)=ε(λ),Tt=Tb=T,則公式(6)簡化為

(8)

只要獲取吸收通道內各波長處的大氣效應參數D(λ)和S(λ),結合地表發射率和溫度數據,即能計算得到各波長處的地表輻亮度LBOA(λ); 然后,利用通道響應函數(圖1)對LBOA(λ)進行通道等效計算,即可模擬得到S1和S2通道的地表輻亮度。

1.2.2 模型大氣效應參數計算方法

對于原始輻射傳輸解析模型,Cota等通過MODTRAN[17]獲取同一大氣條件下2組不同地表反射率的入瞳輻亮度,計算得到這一大氣狀況下的9個大氣效應參數[15]。同理,該方法也適用于D(λ)和S(λ)的獲取。

首先,將MODTRAN輸入參數中的傳感器高度設為0.000 1 m,假設地表海拔為海平面,則得到的輻亮度可近似認為是地表輻亮度。此外,地表溫度設為0.001 K,用于忽略地表發射,即

(9)

此時,計算同一大氣條件下2組不同反射率ρ1和ρ2的輻亮度,得到的地表反射輻亮度結果為

(10)

(11)

計算得到S(λ)和D(λ),即

(12)

(13)

地表輻亮度LBOA1(λ)和LBOA2(λ)從MODTRAN的輸出結果獲得[17]。因此,基于MODTRAN計算的LBOA1(λ)和LBOA2(λ),通過(11)(12)獲得大氣效應參數后,結合地表發射率和溫度信息,以及通道響應函數即可模擬4.3 μm吸收通道的地表輻亮度。

2 模擬實驗和分析

2.1 地表輻亮度圖像模擬

2.1.1 地表發射率圖像模擬

選擇中國東北地區的MODIS23通道的地表數據作為實驗區,該圖像獲取于2004年10月10日,空間分辨率為1 km,實驗區為240像元×240像元。基于MODIS23通道的地表發射率圖像,通過轉換模型公式(2)(3)獲得模擬通道S1和S2的地表發射率圖像,如圖4所示。

(a) S1 (b) S2

圖4SPIRIT-Ⅲ波段S1和S2的地表發射率模擬圖像

Fig.4SimulatedsurfaceemissivityimagesofSPIRIT-ⅢbandS1andS2

2.1.2 地表輻亮度圖像模擬

地表溫度圖像采用對應的MODIS地表溫度產品MOD11_L2。

考慮到MODIS圖像幅寬大,像元間觀測幾何差異明顯,而大氣效應參數受觀測角度影響會發生變化,因此需要逐像元求解D(λ)和S(λ)。鑒于地表輻亮度主要受太陽天頂角的影響[9],本文運用查找表法,預先計算一系列太陽天頂角下的大氣效應參數,然后依據MOD03產品提供的太陽天頂角圖像,基于插值算法獲取每個像元對應觀測角度下的D(λ)和S(λ)。建立D(λ)和S(λ)的查找表時,根據實驗區的數據獲取時間和經緯度,選擇標準亞極地夏季大氣作為建立查找表的大氣模式。查找表建立時的其他主要輸入參數如下: 大氣類型為標準亞極地夏季; 氣溶膠類型為鄉村型; 水汽含量為2.08 g/cm2; 能見度為23 km; 太陽天頂角為0~70°,以5°為間隔。

綜上所述,基于公式(8),結合發射率模擬圖像、地表溫度圖像以及插值獲得的每個像元的大氣效應參數,并通過通道等效計算,模擬獲得4.3 μm吸收通道S1和S2的地表輻亮度圖像,如圖5所示。

(a) S1 (b) S2

圖5SPIRIT-Ⅲ波段的地表輻亮度模擬圖像

Fig.5Simulatedbottom-of-atmosphereradianceimageinSPIRIT-ⅢbandsS1andS2

2.2 模擬方法精度分析

由于SPIRIT-Ⅲ數據具有紅外預警的用途,其4.3 μm強吸收通道的圖像無法獲取到。再者,即使輻亮度圖像能夠獲取,在大氣強吸收通道也難以通過大氣糾正獲得高精度的地表發射率圖像。因此,本文采用以下方式對提出的模擬方法進行精度驗證: ①將MODIS發射率光譜庫中除樣本數據以外的23種地物作為檢驗數據,評價發射率通道轉換模型的模擬精度; ②通過與MODTRAN計算結果比較,評價地表輻亮度計算方法的模擬精度。

2.2.1 發射率轉換模型精度分析

基于公式(2)(3)轉換得到用于模型驗證的23種地物在S1/S2通道的模擬發射率。將模擬發射率與實際發射率比較,計算相對誤差,結果見圖6。

圖6 模擬的通道等效發射率的相對誤差

對于這23種地物,S1通道發射率模擬的相對誤差為-4.28%~5.99%,S2通道發射率模擬的相對誤差為-5.99%~4.78%。由此可見,采用通道轉換模型能夠準確模擬吸收通道的地表發射率。

2.2.2 地表輻亮度模擬精度分析

對于地表輻亮度模擬,分別計算不同的大氣模式和觀測幾何條件下,基于本文的地表輻亮度解析模型計算的模擬地表輻亮度,和MODTRAN直接計算的地表輻亮度。視后者為真值,計算前者模擬地表輻亮度的相對誤差,從而評價本文的地表輻亮度模擬方法能否替代MODTRAN的計算。模擬計算的主要輸入參數: 大氣類型包括熱帶型(TRP)、中緯度夏季型(MLS)、中緯度冬季型(MLW)、亞極地夏季型(SAS)、亞極地冬季型(SAW); 太陽天頂角(Solar zenith angle,SZA)分別為0°, 30°和60°; 觀測天頂角為30°; 相對方位角為90°。相對誤差結果見圖7。

(a) S1(b) S2

圖7S1和S2波段的模擬地表輻亮度的相對誤差

Fig.7Relativeerrorsofsimulatedbottom-of-atmosphereradianceinbandsS1andS2

圖7表明,在5種標準大氣和3種太陽角度條件下,采用解析模型替代大氣輻射傳輸模型MODTRAN計算SPIRIT-Ⅲ兩個通道的地表輻亮度,具有極高的精度,相對誤差不超過±0.02%。由此可見,采用本文提出的方法進行吸收通道的地表輻亮度模擬可行。

3 結論

1)本文發展了一種利用MODIS23通道數據進行以4.3 μm為中心吸收波長的大氣強吸收通道的地表輻亮度圖像模擬方法。該方法主要包括2個步驟: 首先,基于鄰近通道地表發射率之間的強相關性,建立發射率通道轉換模型,模擬吸收通道的地表發射率圖像; 然后,建立4.3 μm吸收通道地表輻亮度的解析模型,用于實現大氣效應參數的計算和地表輻亮度圖像的模擬。

2)由于真實吸收通道數據的獲取困難,本文分別基于光譜庫數據和輻射傳輸模型MODTRAN計算結果,對用于地表發射率模擬的通道轉換模型和地表輻亮度解析模型的模擬結果進行了精度評價。相對誤差顯示地表發射率模擬誤差不超過±6%; 地表輻亮度的模擬誤差則不超過±0.02%。這表明本文提出的地表發射率和地表輻亮度計算模型均能實現吸收通道圖像的精確模擬。

3)本文方法將為4.3 μm吸收通道的入瞳輻亮度模擬、遙感器成像過程模擬等后續工作提供數據支撐。此外,如果需要模擬高空間分辨率的4.3 μm吸收通道圖像時,地表模擬不能忽略地物BRDF特性的影響,這部分內容將在后續研究中進一步加強。

[1] 葉慶,孫曉泉,邵立.紅外預警衛星最佳探測波段分析[J].紅外與激光工程,2010,39(3):389-393. Ye Q,Sun X Q,Shao L.Analysis of optimum detective wavebands for infrared early-warning satellite[J].Infrared and Laser Engineering,2010,39(3):389-393.

[2] 蔣躍,鄧磊,臧鵬.美國天基紅外預警系統的發展現狀和技術特點[J].空軍雷達學院學報,2011,25(2):105-108. Jiang Y,Deng L,Zang P.Developmental state and technical features of American space-based infrared early warning system[J].Journal of Air Force Radar Academy,2011,25(2):105-108.

[3] 陳方,牛錚,廖楚江.遙感圖像模擬技術方法與應用分析[J].地球信息科學,2006,8(3):114-118. Chen F,Niu Z,Liao C J.Analysis on simulation of remote sensing image and its application[J].Geo-Information Science,2006,8(3):114-118.

[4] 孫偉健,林軍,阮寧娟,等.國外光學遙感成像系統仿真軟件發展綜述與思考[J].航天返回與遙感,2010,31(3):70-75. Sun W J,Lin J,Ruan N J,et al.Summarization and consideration of oversea’s simulation software development for optical remote sensing system[J].Spacecraft Recovery and Remote Sensing,2010,31(3):70-75.

[5] 馬曉珊,孟新,楊震,等.光學遙感成像系統全鏈路仿真框架研究[J].量子電子學報,2012,29(4):392-399. Ma X S,Meng X,Yang Z,et al.Framework of entire image chains simulation for optical remote sensing images system[J].Chinese Journal of Quantum Electronics,2012,29(4):392-399.

[6] 徐大琦,杜永明,林軍,等.陸地觀測衛星圖像模擬技術研究[J].國土資源遙感,2015,27(2):80-87.doi:10.6046/gtzyyg.2015.02.13. Xu D Q,Du Y M,Lin J,et al.Research on image simulation technology of land observation satellite[J].Remote Sensing for Land and Resources,2015,27(2):80-87.doi:10.6046/gtzyyg.2015.02.13.

[7] Bartschi B Y,Morse D E,Woolston T L.The spatial infrared imaging telescope Ⅲ[J].Johns Hopkins APL Technical Digest,1996,17(2):215-225.

[8] Stair A T.MSX design parameters driven by targets and backgrounds[J].Johns Hopkins APL Technical Digest,1996,17(1):11-18.

[9] Zhang B,Liu Y,Zhang W J,et al.Analysis of the proportion of surface reflected radiance in mid-infrared absorption bands[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2014,7(6):2639-2646.

[10]Li Z L,Tang B H,Wu H,et al.Satellite-derived land surface temperature:Current status and perspectives[J].Remote Sensing of Environment,2013,131:14-37.

[11]Reda I,Andreas A.Solar position algorithm for solar radiation applications[J].Solar Energy,2004,76(5):577-589.

[12]Demir B,Celebi A,Erturk S.A low-complexity approach for the color display of hyperspectral remote-sensing images using one-bit-transform-based band selection[J].IEEE Transactions on Geoscience and Remote Sensing,2009,47(1):97-105.

[13]Zhang Y L.MODIS UCSB Emissivity Library[EB/OL].(1999-11-10)[2015-07-24].http://www.icess.ucsb.edu/modis/EMIS/html/em.html.

[14]Egan M P,Price S D,Moshir M M,et al.The Midcourse Space Experiment Point Source Catalog Version 1.2 Explanatory Guide[R].Hanscom Afb,MA:Air Force Research Lab,1999.

[15]Cota S A,Kalman L S.Predicting top-of-atmosphere radiance for arbitrary viewing geometries from the visible to thermal infrared:Generalization to arbitrary average scene temperatures[C]//Proceedings SPIE 7813,Remote Sensing System Engineering Ⅲ.San Diego,California,United States:SPIE,2010:781307.

[16]Cota S A,Bell J T,Boucher R H,et al.PICASSO:An end-to-end image simulation tool for space and airborne imaging systems[J].Journal of Applied Remote Sensing,2010,4(1):043535.

[17]Berk A,Anderson G P,Acharya P K,et al.MODTRAN4 User’s Manual[M].Hanscom Afb,MA:Air Force Research Laboratory,Space Vehicles Directorate,Air Force Materiel Command,1999.

(責任編輯:邢宇)

Radianceimagesimulationatthebottomofatmosphereinmid-infraredabsorptionbands

LIU Yao1, ZHANG Wenjuan2, ZHANG Bing2, GAN Fuping1

(1.ChinaAeroGeophysicalSurveyandRemoteSensingCenterforLandandResources,Beijing100083,China; 2.KeyLaboratoryofDigitalEarthScience,InstituteofRemoteSensingandDigitalEarth,ChineseAcademyofSciences,Beijing100094,China)

This method was illustrated by applying it to simulating bottom-of-atmosphere(BOA) radiance for two 4.3 μm absorption bands of the SPIRIT-Ⅲ sensor by using MODIS data of band 23, which is close to 4.3 μm. First, surface emissivity images in these two 4.3 μm absorption bands were simulated using band translation models. Second, analytic model of BOA radiance was deduced based on an existing analytic model in mid-infrared bands, and then it was combined with simulations from radiative transfer model MODTRAN to calculate parameters of the atmospheric effects for these 4.3 μm absorption bands. Finally, based on the proposed analytic model, BOA radiance in SPIRIT-Ⅲ’s two absorption bands can be generated from surface emissivity, temperature, atmospheric effect parameters and SPIRIT-Ⅲ’s spectral response functions. Accuracy assessment on the simulation results shows that this method can produce surface emissivity and BOA radiance with errors less than 6% and 0.02%, respectively. Therefore, the method proposed in this paper can effectively and precisely simulate BOA radiance for the 4.3 μm absorption bands, and provide radiance datasets for the at-sensor radiance simulation and sensor imaging simulation.

image simulation; mid-infrared; atmospheric absorption; surface emissivity

10.6046/gtzyyg.2017.03.14

劉瑤,張文娟,張兵,等.中紅外大氣強吸收通道的地表輻亮度圖像模擬[J].國土資源遙感,2017,29(3):98-103.(Liu Y,Zhang W J,Zhang B,et al.Radiance image simulation at the bottom of atmosphere in mid-infrared absorption bands[J].Remote Sensing for Land and Resources,2017,29(3):98-103.)

2015-07-24;

2015-11-20

國家自然科學基金“高光譜遙感智能觀測機理與信息處理模型研究”(編號: 41325004)、國家自然科學基金“高光譜遙感在軌替代定標模型與方法研究”(編號: 41271370)、國家科技重大專項“高分辨率對地觀測系統”(編號: 30-Y20A06-9003-15/16)和中國地質調查局地質調查項目“全國礦產資源開發環境遙感監測(編號: 12120120300016009)共同資助。

劉瑤(1988-),女,博士,工程師,主要從事遙感圖像仿真模擬研究。Email: yao.liu_agrs@foxmail.com。

張文娟(1983-),女,博士,高級工程師,主要從事地表輻射傳輸建模和衛星圖像模擬等研究。Email: zhangwj@radi.ac.cn。

TP 751.1

: A

: 1001-070X(2017)03-0098-06

猜你喜歡
大氣模型
一半模型
大氣的呵護
軍事文摘(2023年10期)2023-06-09 09:15:06
太赫茲大氣臨邊探測儀遙感中高層大氣風仿真
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
大氣古樸揮灑自如
大氣、水之后,土十條來了
新農業(2016年18期)2016-08-16 03:28:27
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
莊嚴大氣 雄媚兼備
主站蜘蛛池模板: 欧美亚洲第一页| 97超爽成人免费视频在线播放| 欧美日韩另类在线| 国产激情在线视频| 国产一在线观看| 五月婷婷欧美| 亚洲天堂777| 视频一区视频二区中文精品| 欧美激情网址| 午夜国产精品视频黄 | 欧美中日韩在线| 国产91在线|日本| 男女男精品视频| 亚洲高清日韩heyzo| 亚洲五月激情网| 玖玖精品视频在线观看| 综合五月天网| 性欧美精品xxxx| 久久精品无码一区二区日韩免费| 亚洲日韩国产精品无码专区| 青青草原国产精品啪啪视频| 熟妇人妻无乱码中文字幕真矢织江| av大片在线无码免费| 国产亚洲精品无码专| 99热这里只有精品国产99| 国产黑丝一区| 亚洲AV无码久久天堂| 国产91透明丝袜美腿在线| 国产青榴视频在线观看网站| 国产麻豆aⅴ精品无码| 任我操在线视频| 18禁黄无遮挡网站| 欧美性色综合网| 久久精品一卡日本电影| 中文一区二区视频| 丁香六月激情综合| 青草视频在线观看国产| 亚洲人在线| 久久综合成人| 国产AV毛片| 91视频国产高清| 免费看美女自慰的网站| 丁香五月亚洲综合在线| 国产午夜不卡| 国产成人亚洲无吗淙合青草| 亚洲综合极品香蕉久久网| 国产精品分类视频分类一区| 不卡网亚洲无码| 精品人妻一区无码视频| 欧美区国产区| a级高清毛片| 国产男女免费视频| 国产一级在线观看www色| 国产成人你懂的在线观看| 99久久国产精品无码| 国产欧美精品午夜在线播放| 午夜精品福利影院| 国产在线97| 污视频日本| 97视频在线观看免费视频| 亚洲VA中文字幕| 中文字幕在线欧美| 欧美日韩国产综合视频在线观看| A级毛片高清免费视频就| 国产在线一二三区| 亚洲成人黄色在线| 美女扒开下面流白浆在线试听 | 国产99精品视频| 欧美日韩v| 国产成人精品2021欧美日韩 | 日韩免费毛片视频| 91在线精品麻豆欧美在线| 国产AV无码专区亚洲精品网站| 国产第一页第二页| 日韩欧美国产区| 波多野结衣一区二区三区AV| 国产亚洲视频免费播放| 澳门av无码| 国产凹凸一区在线观看视频| 欧美日韩成人在线观看| 无码人中文字幕| 国产日产欧美精品|