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

基于ETM+數據的煤田火區溫度異常信息提取

2017-04-27 09:31:56張春森徐肖雷陳越峰
自然資源遙感 2017年2期

張春森, 徐肖雷, 陳越峰

(1.西安科技大學測繪科學與技術學院,西安 710054; 2.新疆煤田滅火工程局,烏魯木齊 830063)

基于ETM+數據的煤田火區溫度異常信息提取

張春森1, 徐肖雷1, 陳越峰2

(1.西安科技大學測繪科學與技術學院,西安 710054; 2.新疆煤田滅火工程局,烏魯木齊 830063)

快速而準確地提取煤田火區溫度異常信息對于主動發現煤火和及時治理煤火具有重要意義。以新疆維吾爾自治區烏魯木齊市郊的大泉湖火區為研究區,基于ETM+遙感影像數據,采用普適性單通道法反演煤田火區地表溫度,利用人工閾值法和密度分割法提取背景區與溫度異常區,將溫度異常區與磁法成果圖進行圖層疊加分析。研究結果表明,普適性單通道算法反演的溫度均方根誤差為0.68℃,溫度異常區與煤火區的重疊率為82.71%,溫度異常區反演準確率為80.17%,為基本圈定煤火范圍提供了有力的參考。

煤田火區; 溫度反演; 溫度異常; 火區圈定; 磁法勘探

0 引言

煤田火災是一個世界性的自然災害,不僅會造成巨大的經濟損失,還會造成嚴重的環境污染。遙感技術具有大范圍、全天候、實時、快速的特點,利用遙感技術來探測煤火和初步勘察火區狀況是煤炭災害領域的一個研究方向,也是煤田滅火工程開展的基礎。利用遙感方法進行監測煤火方面,國內外學者已做了大量的研究工作,早在1964年Slavecki[1]就開始利用熱紅外相機探測地下煤火,國內學者覃志豪等[2]針對僅有一個熱紅外波段的Landsat TM數據提出了單窗算法; 毛克彪等[3]、張敦虎等[4]、張曉等[5]先后進行了劈窗算法的研究; 蔣大林等[6]、徐涵秋等[7]、王亞維等[8]分別利用Landsat 8 TIRS和FY-2C數據進行了地表溫度的反演試驗; 張秀山[9]、蔡忠勇等[10]對新疆煤田火災的現狀進行了分析研究; 張躍彬[11]、劉碩等[12]對煤火范圍的圈定進行了研究。雖然相關研究已做了很多,但仍存在一些問題,如利用Landsat 8 TIRS數據進行溫度反演受大氣水汽含量的影響較大,且空間分辨率低,僅為100 m,局限性明顯。煤火范圍的圈定僅考慮重疊率一項評價指標,忽視了方法自身的準確率評價等。

本文利用ETM+數據,采用普適性單通道算法反演煤火區地表溫度,利用人工閾值法和密度分割法提取溫度異常區,將溫度異常區與磁法成果圖進行疊加分析,計算二者的重疊率和反演的準確率,為采用遙感方法探測煤火和初步圈定煤火范圍提供技術支撐。

1 研究方法

1.1 煤火區地表溫度反演

單通道算法由Jimenez等[13]提出的基于一個熱紅外通道利用遙感影像反演地表溫度的適普性算法,其計算公式為

T=γ[ε-1(ψ1Lsensor+ψ2)+ψ3]+δ,

(1)

式中:T為反演得到的地表溫度;ε為地表比輻射率;Lsensor為星上輻射亮度;γ和δ為Planck函數進行泰勒級數展開時的相關參數,其計算公式為:

(2)

δ=-γLsensor+Tsensor,

(3)

Tsensor=K2/ln(K1/Lsensor+1),

(4)

Lsensor=(Lmax-Lmin)/255DN+Lmin。

(5)

式中:Tsensor為星上輻射亮度對應的亮度溫度;c1和c2為普朗克常量,c1=1.19 104×108W·μm4·m-2·sr-1,c2=1.43 877×104μm·K;λ為熱紅外波段的有效波長;K1和K2為熱紅外波段的定標常數,K1=666.09 W·m-2·sr-1·μm-1,K2=1 282.71K;Lmax和Lmin分別為熱紅外波段光譜輻射值的上限和下限;DN為像元亮度值。

式1中ψ1,ψ2,ψ3為大氣函數,其公式為

(6)

式中:pij(i,j=1,2,3)為與水汽含量ω相關的大氣參數;ω為大氣水汽含量,由飽和水汽壓、大氣溫度和濕度計算得到。

在計算地表比輻射率之前,需首先對研究區進行監督分類,根據研究區的地貌特點,將影像分為植被、裸地及植被裸地混合區3種地物類型,并計算其歸一化植被指數(nermalized difference vegetation index,NDVI)值。具體計算公式為

Fv=[(NDVI-NDVIs)/(NDVIv-NDVIs)]2,

(7)

(8)

式中:Fv為植被覆蓋度;NDVI為歸一化植被指數,NDVI和NDVIv分別為裸地區和植被區的NDVI值;εv和εs分別為植被區和裸地區比輻射率;ρ為紅光波段的反射率;dε為地表比輻射率改正值。

當NDVINDVIv時,視為植被,按照公式(8)中第3式計算,可得εv; 當NDVIs≤NDVI≤NDVIv時,視為植被裸地混合區,且εs和εv已知,按照公式(8)中第2式計算,可以得到混合區的比輻射率。dε是由地表的起伏形態所引起的比輻射率的改正值,由于該實驗區地形較為平坦,故可忽略不計[14]。

1.2 溫度異常區提取

煤層在燃燒過程中,熱量以熱輻射的方式,通過煤火上層巖石的熱傳導和地層裂隙中空氣的熱對流向上傳導,在地面和低空中形成溫度相對高于其周圍環境溫度的區域,即為溫度異常區[15]。溫度異常區的提取是探測煤火和圈定煤火范圍的前提。

通過遙感反演得到的煤田火區地表溫度通常受很多因素的影響,如影像成像的時間(白天或夜間)、季節、氣候、燃燒情形等條件的不同都會影響反演結果。如果設置固定的溫度閾值來提取不同時期、不同環境下的溫度異常區,就會產生比較大的差異,不能真實反映煤田火區的狀況。因此有必要針對不同情況下的溫度反演結果,按照統一的方法來計算溫度閾值。首先采用密度分割法將煤火區溫度劃分為不同的溫度區間,然后采用人工閾值法確定提取溫度異常區的閾值,進而區分出溫度異常區與背景區。蔣衛國[16]等研究發現,地表溫度均值與2倍標準偏差之和是溫度異常區與背景區的最佳閾值分割點。若溫度高于設定的閾值,則為溫度異常區; 反之,則為背景區。計算公式為

(9)

式中:Ti為地表溫度圖像中任一像元的地表溫度值;N為地表溫度圖像的像元總數;Ta為地表溫度的平均值;Tδ為地表溫度標準偏差;T閾為分割點閾值。

1.3 溫度異常區面積估算

在提取到溫度異常區之后,利用ENVI軟件統計分析模塊中的快速統計(quick statistics)功能,對溫度異常區影像進行分析。該功能可以顯示影像任一像元的溫度,溫度出現的頻率等信息。通過像元與溫度的對應關系,結合溫度出現的頻率和影像的空間分辨率,可以快速準確地計算出溫度異常區的面積,即

(10)

式中:Sr為溫度異常區面積,m2;nTi為溫度閾值范圍內某一溫度Ti對應的像元個數;u為溫度閾值范圍內不同溫度的個數;Q為ETM+影像B61波段的空間分辨率(60 m)。

1.4 技術路線與算法實現

地表溫度反演之前,需對ETM+影像進行輻射定標和大氣校正等預處理。輻射定標是計算星上輻射亮度Lsensor的過程,即利用ETM+影像各波段的增益值和偏置值計算相應波段的輻射亮度。大氣校正是利用大氣校正工具來消除大氣的影響。然后結合1∶2 000地形圖對ETM+影像進行幾何精糾正,最后將校正好的影像作為地表溫度反演的輸入數據。

根據煤田火區地表溫度反演的算法及流程,利用IDL數組操作函數進行波段運算,實現地表溫度的反演和溫度異常區的提取。技術流程如圖1所示。

圖1 技術流程圖

利用影像DN值和ETM+定標系數計算星上輻射亮度Lsensor,利用星上輻射亮度計算亮度溫度Tsensor,利用歸一化植被指數和植被覆蓋度Fv計算地表比輻射率ε; 再利用中間變量γ,δ和ψ等參數計算煤火區地表溫度T,并統計地表溫度平均值Ta和標準偏差Tδ,計算分割點閾值,提取溫度異常區; 然后將溫度異常區與磁法成果圖進行疊加分析,最后初步圈出煤火分布范圍。

2 數據源及實驗分析

2.1 實驗區及數據源概況

本文選擇大泉湖煤火區為實驗區。該區位于新疆維吾爾自治區烏魯木齊市境內,海拔1 006 m,屬溫帶大陸性干旱氣候,地理坐標為E 87°25′29.71″~87°27′17.25″, N 43°47′31.23″ ~43°47′57.15″,面積為2.5 km2左右,是新疆煤田火區中比較典型的代表。

實驗選擇2009年8月25日凌晨獲取的1景ETM+第6波段熱紅外數據(波長范圍為10.40~12.50 μm,空間分辨率為60 m),用于反演地表溫度,提取溫度異常區。當天的氣象數據來自新疆烏魯木齊氣象局,包括大氣溫度、濕度、飽和水汽壓等數據,用來計算大氣水汽含量和進行大氣校正。新疆煤田滅火工程局實測的2009年8月25日地表溫度數據用來對比驗證ETM+反演的地表溫度。2007年5月的煤火區磁法勘探數據用來與溫度異常區進行疊加分析。2011年2月的煤火區地1∶2 000形圖用來分析地形和進行幾何精糾正。

2.2 溫度反演結果與分析

根據大泉湖火區的ETM+熱紅外遙感數據,按照上述算法及流程,反演得到了大泉湖火區的地表溫度分布(圖2)。其中反演得出的地表最低溫度為23.3 ℃,最高為57.2 ℃。

圖2 反演地表溫度分布圖

為驗證反演結果,通常采取與實測衛星過境時的地表溫度或利用大氣模型模擬得到的地表溫度進行對比的方法。本實驗采用由新疆煤田滅火工程局提供的衛星過境時的地面實測數據。衛星過境時間為地方時10時28分。在衛星過境的前后2 h內,進行實驗區地表溫度的采集。利用手持式ZyTemp TN435紅外測溫儀和GPS-RTK進行地表溫度的網格式測量,即2個人分別從實驗區的同側2個角點出發,相向而行,以20 m為間隔,之字線路依次觀測,記錄坐標值與溫度值。取60 m×60 m大小網格(與ETM+熱紅外數據的地面分辨率相同)中溫度平均值作為實測溫度值。均勻選取10個實測溫度值,與反演的相應地表溫度進行回歸分析[17],結果如圖3所示。

圖3 反演溫度與實測溫度回歸分析圖

圖3表明普適性單通道算法反演的地表溫度值與實測地表溫度值成一元線性相關,其相關系數為0.955 9,相關性較好。反演得到的地表溫度值與實測地表溫度值間的均方根誤差為0.681 2 ℃。

2.3 溫度異常區提取結果與分析

采用人工閾值法和密度分割法提取研究區的溫度異常區。由公式(9)計算的溫度閾值為53.8 ℃,因此將背景區溫度設定為23.3 ℃~53.8 ℃,溫度異常區溫度設定為53.8 ℃~57.2 ℃,再次對地表溫度分布圖進行密度分割。按照上述算法及流程,得到大泉湖火區的溫度異常分布,如圖4所示。

圖4 實驗區溫度異常區分布圖

為驗證溫度異常區地理位置分布的準確性,采用新疆煤田滅火工程局物探部門提供的大泉湖火區磁法勘探數據進行驗證。該部門的磁法勘探按照工作規范要求在E 87°25′07.43″ ~ 87°27′17.23″,N 43°47′30.92″~ 43°47′57.19″范圍內布設了測線,根據探測對象與圍巖之間的磁性差異,利用質子磁力儀觀測記錄磁異常,通過外業測量和內業解算,繪得火區的燃燒邊界與燃燒深度。

多年來,許多學者的研究已經證明,磁法勘探可以有效地確定著火點位置和圈定火區范圍, 并利用磁異常的變化情況大致判斷煤火走勢[11]。因此磁測方法目前已成為煤田滅火工程部門地面煤火探測的有效方法之一。為此,本文選取實驗區內磁測范圍作為重點區域進行對比驗證(圖5所示)。

圖5 磁法探測區位置

將反演提取的溫度異常與磁測成果圖通過圖層疊加,如圖6所示。計算二者的重疊率和溫度異常區反演的準確率,其中:

(11)

圖6 圖層疊加效果

(12)

式中:RO為重疊率;Ra為準確率;Sm為磁法所測的燃燒面積(視為相對真值),即燃燒邊界圈定的范圍(圖7所示);Sr為遙感反演的溫度異常區面積;SO為Sm與Sr重疊部分的面積。具體數值見表1。

圖7 圖層疊加局部放大

表1 面積對比一覽表

由表1可知,計算的溫度異常區面積Sr為0.417 6 km2,落入磁測燃燒邊界內的溫度異常區(重疊部分)面積SO為0.334 8 km2; 磁法所測的燃燒范圍Sm為0.404 8 km2。由公式(11)、(12)計算得溫度異常區與磁法所測燃燒范圍的的重疊率為82.71%,溫度異常區反演的準確率為80.17%。

基于遙感影像地表溫度反演提取的溫度異常區與磁法勘探確定的燃燒范圍出現差別的原因,可歸納為2個方面: 一是磁法勘探數據的時間為2007年,而ETM+遙感數據的時間為2009年,在此期間,火區內隨著煤火的蔓延和擴展,出現了一些新的著火煤層露頭,導致反演得到的部分溫度異常區出現在磁測燃燒邊界之外; 相反,火區內部分老的著火點可能已經熄滅,導致磁測燃燒邊界內反演結果出現背景區。二是ETM+熱紅外波段的空間分辨率有限,導致部分零星溫度異常點沒有被反演出來而歸為背景區。

3 結論

1)普適性的熱紅外單通道地表溫度反演算法應用于煤田火區的地表溫度反演是可行的,且精度較高,溫度誤差可控制在1℃以內。該方法對于主動發現自燃煤田的溫度異常區,及時治理煤火具有重要意義。

2)人工閾值法和密度分割法可以準確高效地提取地表溫度異常區,且提取效果較好,重疊率和準確率均在80%以上,該方法對于初步圈定煤火范圍,監測煤火火情變化具有指導意義。

通過遙感地質解譯方法進行煤田火區構造裂縫和導火帶的提取與分析,實現火區的綜合監測與防治將是我們今后的研究工作。

致謝: 本研究得到了新疆煤田滅火工程局的大力支持和幫助,感謝新疆煤田滅火工程局陳越峰、曾強等提供的地表溫度實測數據與磁法勘探數據,再次表示感謝!

[1] Slaveck R J.Detection and location of subsurface coal fires[C]//Proceedings of the 3rd International Symposium on Remote Sensing of Environment.Ann Arobor,Michigan,U.S.A.:University of Michigan,1964:537-547.

[2] 覃志豪,張明華,Karnieli A,等.用陸地衛星TM6數據演算地表溫度的單窗算法[J].地理學報,2001,56(4):456-466. Qin Z H,Zhang M H,Karnieli A,et al.Mono-window algorithm for retrieving land surface temperature from Landsat TM6 data[J].Acta Geographica Sinica,2001,56(4):456-466.

[3] 毛克彪,施建成,覃志豪,等.一個針對ASTER數據同時反演地表溫度和比輻射率的四通道算法[J].遙感學報,2006,10(4):593-599. Mao K B,Shi J C,Qin Z H,et al.A four-channel algorithm for retrieving land surface temperature and emissivity from ASTER data[J].Journal of Remote Sensing,2006,10(4):593-599.

[4] 張敦虎,王曉鵬,康高峰,等.基于ASTER數據的煤田火區溫度異常信息提取方法研究——以新疆輪臺縣陽霞火區為例[J].中國煤炭地質,2009,21(9):35-37. Zhang D H,Wang X P,Kang G F,et al.Study on coal fire area land surface temperature anomaly information extraction from ASTER satellite images:A case study of Yangxia coal fire area in Luntai County,Xinjiang[J].Coal Geology of China,2009,21(9):35-37.

[5] 張曉,湯瑜瑜,黃小仙,等.用8.0~9.3μm遙感數據反演地溫的分裂窗算法性能分析[J].國土資源遙感,2015,27(2):88-93.doi:10.6046/gtzyyg.2015.02.14. Zhang X,Tang Y Y,Huang X X,et al.Performance analysis of split-window algorithms for retrieving land surface temperature using remote sensing data of 8.0~9.3μm[J].Remote Sensing for Land and Resources,2015,27(2):88-93.doi:10.6046/gtzyyg.2015.02.14.

[6] 蔣大林,匡鴻海,曹曉峰,等.基于Landsat 8的地表溫度反演算法研究——以滇池流域為例[J].遙感技術與應用,2015,30(3):448-454. Jiang D L,Kuang H H,Cao X F,et al.Study of land surface temperature retrieval based on Landsat 8-with the sample of Dianchi Lake Basin[J].Remote Sensing Technology and Application,2015,30(3):448-454.

[7] 徐涵秋,林中立,潘衛華.單通道算法地表溫度反演的若干問題討論——以Landsat系列數據為例[J].武漢大學學報:信息科學版,2015,40(4):487-492. Xu H Q,Lin Z L,Pan W H.Some issues in land surface temperature retrieval of Landsat thermal data with the single-channel algorithm[J].Geomatics and Information Science of Wuhan University,2015,40(4):487-492.

[8] 王亞維,宋小寧,唐伯惠,等.基于FY-2C數據的地表溫度反演驗證——以黃河源區瑪曲為例[J].國土資源遙感,2015,27(4):68-72.doi:10.6046/gtzyyg.2015.04.11. Wang Y W,Song X N,Tang B H,et al.Validation of FY-2C derived land surface temperature over the source region of the Yellow River:A case study of Maqu country[J].Remote Sensing for Land and Resources,2015,27(4):68-72.doi:10.6046/gtzyyg.2015.04.11.

[9] 張秀山.新疆煤田火燒區特征及滅火問題探討[J].中國煤田地質,2004,16(1):18-21. Zhang X S.Probe into Xinjiang coalfield underground combustion area features and fire-fighting problems[J].Coal Geology of China,2004,16(1):18-21.

[10]蔡忠勇,魏軍.新疆煤田火災現狀及其應對措施[J].煤礦安全,2008,34(10):93-94. Cai Z Y,Wei J.Coal field fires in Xinjiang:Its current state, challenges and countermeasures[J].Safety in Coal Mines,2008,34(10):93-94.

[11]張躍彬.磁法勘探在圈定煤田自燃區及在煤田地質勘察中的應用研究[J].中國新技術新產品,2009(3):93. Zhang Y B.Research of magnetic method applicated in coalfield geological exploration and delineating coalfield area[J].China New Technologies and Products,2009(3):93.

[12]劉碩,王涌宇,王俊峰,等.基于地表異常特征的淺埋藏煤層火區邊界劃定[J].煤炭科學技術,2014,42(S1):132-133,136. Liu S,Wang Y Y,Wang J F,et al.Fire area boundaries determination of shallow depth coal seam based on surface abnormal features[J].Coal Science and Technology,2014,42(S1):132-133,136.

[14]Sobrino J A,Jimenez-Muoz J C,Soria G,et al.Land surface emissivity retrieval from different VNIR and TIR sensors[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(2):316-327.

[15]姬洪亮.多尺度熱紅外遙感數據在煤田火區信息提取中的應用[D].烏魯木齊:新疆大學,2012. Ji H L.Application of Multi-Scale Thermal Infrared Remote Sensing Data in Coalfield Fire Area Information Extraction[D].Urumqi:Xinjiang University,2012.

[16]蔣衛國,武建軍,顧磊,等.基于夜間熱紅外光譜的地下煤火監測方法研究[J].光譜學與光譜分析,2011,31(2):357-361. Jiang W G,Wu J J,Gu L,et al.Monitoring method of underground coal fire based on night thermal infrared remote sensing technology[J].Spectroscopy and Spectral Analysis,2011,31(2):357-361.

[17]姬洪亮,塔西甫拉提·特依拜,蔡忠勇,等.基于TM數據的地下煤火區地表溫度反演與驗證[J].國土資源遙感,2012,24(4):101-106.doi:10.6046/gtzyyg.2012.04.17. Ji H L,Tashpolat·Tiyip,Cai Z Y,et al.The inversion and verification of land surface temperature for coal fire areas based on TM data[J].Remote Sensing for Land and Resources,2012,24(4):101-106. doi:10.6046/gtzyyg.2012.04.17.

(責任編輯: 李瑜)

Temperature anomaly information extraction in coalfield fire area based on ETM+ data

ZHANG Chunsen1, XU Xiaolei1, CHEN Yuefeng2

(1.CollegeofGeomatics,Xi’anUniversityofScienceandTechnology,Xi’an710054,China; 2.XinjiangCoalfieldFireExtinguishingBureau,Urumqi830063,China)

The initiative finding of coal fire and timely governing of coal fire so as to quickly and accurately extract temperature anomaly information are of great importance. In this study, with the coal fire area of Daquanhu in the suburbs of Urumqi of Xinjiang as the study area and ETM+ remote sensing data as the base, the authors used generalized single-channel method to retrieve land surface temperature of coalfield fire area, and then used manual threshold method and density slicing method to extract background area and temperature anomaly. Finally, temperature anomaly area image was superimposed upon the magnetic prospecting resultant map to perform analysis. The results show that the retrieved RMSE of Generalized Single-Channel Method is 0.68℃, the overlap rate between temperature anomaly area and definitized fire area range is 82.71%, and the accuracy rate of temperature anomaly area retrieval is 80.17%. This method can preliminarily delineate the coal fire range and also provides a reference for precisely measuring the coal fire range.

coalfield fire area; temperature retrieve; temperature anomaly; delineate coal fire range; magnetic prospecting

10.6046/gtzyyg.2017.02.29

張春森,徐肖雷,陳越峰.基于ETM+數據的煤田火區溫度異常信息提取[J].國土資源遙感,2017,29(2):201-206.(Zhang C S,Xu X L,Cheng Y F.Temperature anomaly information extraction in coalfield fire area based on ETM+ data[J].Remote Sensing for Land and Resources,2017,29(2):201-206.)

2015-10-16;

2015-12-09

國家自然科學基金“基于圖論的面向對象高分辨率遙感影像多尺度分割方法研究”(編號: 41101410)及陜西省自然科學基金“基于尺度下降技術遙感圖像模擬方法研究”(編號: 2010JM5009)共同資助。

張春森(1963-),男,博士,教授,主要從事攝影測量與遙感應用方面的研究。Email:zhchunsen@aliyun.com。

TP 79

A

1001-070X(2017)02-0201-06

主站蜘蛛池模板: 免费无码AV片在线观看国产| 久久久久九九精品影院| 午夜不卡福利| a欧美在线| 囯产av无码片毛片一级| 日韩麻豆小视频| 国产精品制服| 综合色在线| 日本人真淫视频一区二区三区| 久久精品中文字幕少妇| 国产亚洲美日韩AV中文字幕无码成人| 免费一级无码在线网站| 亚洲国产综合自在线另类| 欧美a在线视频| 国产成人一区在线播放| 国产三级视频网站| 综合色88| 996免费视频国产在线播放| 国产视频一区二区在线观看| 亚洲精品制服丝袜二区| 国产欧美精品专区一区二区| 永久在线精品免费视频观看| 国产综合另类小说色区色噜噜| 超薄丝袜足j国产在线视频| 免费看黄片一区二区三区| 国产一级片网址| 亚洲中字无码AV电影在线观看| 精品久久香蕉国产线看观看gif| 亚洲成人精品在线| 婷婷中文在线| 色老头综合网| 2021精品国产自在现线看| 麻豆精品在线| 国产午夜人做人免费视频中文| 77777亚洲午夜久久多人| 亚洲视屏在线观看| 久久精品只有这里有| 中文字幕自拍偷拍| 高清不卡一区二区三区香蕉| 狠狠色成人综合首页| 又黄又湿又爽的视频| 国产精品亚洲精品爽爽| 国产凹凸视频在线观看| 欧美日韩在线观看一区二区三区| yjizz视频最新网站在线| 欧美日韩另类在线| 五月综合色婷婷| 免费观看男人免费桶女人视频| 国产精品lululu在线观看 | 久久综合亚洲鲁鲁九月天| 免费国产福利| 欧美国产日本高清不卡| 欧美成人综合视频| 久久久久免费精品国产| 亚洲天堂免费在线视频| 毛片免费高清免费| 天堂在线亚洲| 又黄又爽视频好爽视频| 久热中文字幕在线观看| 亚洲一区网站| 国产精品人成在线播放| 亚洲无码91视频| 亚洲第一成年网| 欧美精品v| 1024你懂的国产精品| 99er精品视频| 久久国产V一级毛多内射| 婷婷色一二三区波多野衣| 一级毛片免费不卡在线视频| 中文精品久久久久国产网址| 97久久免费视频| 99一级毛片| 欧美色图久久| 国产理论一区| 国产精品人人做人人爽人人添| 丁香婷婷在线视频| av在线人妻熟妇| 乱色熟女综合一区二区| 夜夜拍夜夜爽| 日韩在线成年视频人网站观看| 亚洲精品动漫| 亚洲人成人无码www|