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

時(shí)序InSAR技術(shù)在大型滑坡監(jiān)測(cè)中的應(yīng)用

2019-01-26 10:25:04趙寶強(qiáng)韓守富白艷萍馬金輝
科技創(chuàng)新與應(yīng)用 2019年1期
關(guān)鍵詞:變形

趙寶強(qiáng) 韓守富 白艷萍 馬金輝

摘 要:綜合PS-InSAR和SBAS-InSAR兩種時(shí)序InSAR技術(shù),對(duì)青海省樂(lè)都區(qū)高家灣滑坡地表變形過(guò)程進(jìn)行遙感監(jiān)測(cè)與分析。通過(guò)SBAS-InSAR技術(shù)獲取該滑坡2003-2010年地表變形數(shù)據(jù),分析歷史變化情況。利用PS-time工具分析高家灣滑坡2014-2017年P(guān)S-InSAR獲得的地表變形監(jiān)測(cè)數(shù)據(jù),將滑坡地表形變過(guò)程劃分為線性和非線性,重點(diǎn)分析了地表形變速率、累積位移突變階段和主裂縫形變情況。研究表明:2003-2010年期間該滑坡地表活動(dòng)較為穩(wěn)定,2014-2017年期間該滑坡發(fā)生了明顯的地表沉降變形且變化趨勢(shì)較為強(qiáng)烈,與滑坡區(qū)實(shí)地勘測(cè)結(jié)果情況一致,較好的證明時(shí)序InSAR技術(shù)用于滑坡監(jiān)測(cè)的有效性。

關(guān)鍵詞:PS-InSAR;SBAS-InSAR;Sentinel-1A;滑坡監(jiān)測(cè)

中圖分類(lèi)號(hào):P237 文獻(xiàn)標(biāo)志碼:A 文章編號(hào):2095-2945(2019)01-0021-04

Abstract: The surface deformation process of Gaojiawan landslide in Ledu District of Qinghai Province's Haidong City is monitored and analyzed using remote sensing in combination with PS-InSAR and SBAS-InSAR time series InSAR technologies. The surface deformation data of the landslide from 2003 to 2010 are obtained by SBAS-InSAR technology, and the historical changes are analyzed. The monitoring data of surface deformation obtained by PS-InSAR in 2014-2017 of Gaojiawan landslide are analyzed using PS-time tool. The process of surface deformation of Gaojiawan landslide is divided into linear and nonlinear ones. The surface deformation rate, the abrupt change stage of cumulative displacement and the deformation of main cracks are emphatically analyzed. The study shows that the surface activity of the landslide is stable in the period 2003-2010, and the landslide has obvious surface subsidence deformation and strong change trend in 2014-2017, which is consistent with the results of the field survey of the landslide area. It is proved that the temporal InSAR technique is effective in landslide monitoring.

Keywords: PS-InSAR; SBAS-InSAR; Sentinel-1A; landslide monitoring

引言

部分斜坡沿著斜坡內(nèi)的一個(gè)或數(shù)個(gè)面在重力作用下做剪切運(yùn)動(dòng)的現(xiàn)象,稱(chēng)為滑坡[1]。滑坡作為一種常見(jiàn)的地質(zhì)災(zāi)害,分布廣,范圍大。我國(guó)西北、西南、華北地區(qū)以及丘陵、黃土高原地區(qū),滑坡的種類(lèi)多、分布面積廣、危害嚴(yán)重且破壞力大,每年都造成巨大的經(jīng)濟(jì)損失,嚴(yán)重制約著災(zāi)害多發(fā)地區(qū)的國(guó)民經(jīng)濟(jì)發(fā)展,威脅著人民生命財(cái)產(chǎn)安全[2]。國(guó)內(nèi)外在監(jiān)測(cè)滑坡災(zāi)害方面的主要手段有:精密水準(zhǔn)測(cè)量、導(dǎo)線測(cè)量、光纖傳感器、全球定位系統(tǒng)GPS、衛(wèi)星遙感技術(shù)、近景攝影測(cè)量、地面激光掃描等[3]。星載合成孔徑雷達(dá)干涉測(cè)量技術(shù)憑借其覆蓋范圍廣、分辨率高、全天時(shí)全天候、監(jiān)測(cè)精度高等特點(diǎn),被廣泛應(yīng)用于地震、火山、冰川移動(dòng)、地下水抽取和地下采礦引起的地表形變監(jiān)測(cè)研究[4]。InSAR(Interferometric Synthetic Aperture Radar)技術(shù)能夠提取地表發(fā)生的微小形變,精度達(dá)到厘米級(jí)甚至是毫米級(jí),其獲取的數(shù)據(jù)具有高精度、高分辨率、覆蓋范圍廣等特點(diǎn)[5]。然而,傳統(tǒng)D-InSAR形變測(cè)量的精度和可靠性受到空間失相關(guān)[6-7]和大氣延遲的嚴(yán)重影響[8-9]。本文結(jié)合地質(zhì)背景資料和野外實(shí)地滑坡勘測(cè)結(jié)果,綜合運(yùn)用PS-InSAR和SBAS-InSAR技術(shù),探討分析時(shí)序InSAR技術(shù)在滑坡變形監(jiān)測(cè)中的適用性。

1 研究區(qū)概況

高家灣滑坡位于青海省海東市樂(lè)都區(qū)洪水鎮(zhèn)湟水河南岸低山丘陵區(qū),為第四系黃土和古近系泥巖組成的復(fù)合型古滑坡。滑坡長(zhǎng)1830m、寬1300m、平均厚110m,規(guī)模2.62×108m3,屬巨型滑坡。2011年至2012年滑坡后壁及后緣發(fā)育長(zhǎng)約500m拉張裂縫。據(jù)2016年1月調(diào)查資料顯示,該裂縫東西延伸增加到950m,其滑體前緣北側(cè)湟水河Ⅱ級(jí)階地表部新增鼓脹裂縫,并使高家灣村局部民房出現(xiàn)變形,嚴(yán)重威脅高家灣村村民及滑坡體中后緣西寧-官亭750千伏輸電線、龍羊峽-海石灣330千伏輸電線、蘭新高速鐵路張家莊隧道及滑坡前緣澀寧蘭輸氣管線、蘭西拉光纜及G109國(guó)道等基礎(chǔ)設(shè)施安全。研究區(qū)中心經(jīng)緯度36°26′35″N,102°32′20″E,南北長(zhǎng)5.24km,東西寬6.12km。地理位置如圖1示。

2 數(shù)據(jù)說(shuō)明及技術(shù)方法

本研究根據(jù)SAR數(shù)據(jù)的分辨率、波長(zhǎng)、觀測(cè)角度、性?xún)r(jià)比和可獲得性等,針對(duì)研究區(qū)地形、植被和變形特征進(jìn)行SAR數(shù)據(jù)源和不同InSAR計(jì)算策略的選擇。Sentinel-1A數(shù)據(jù)周期較短且數(shù)據(jù)豐富,適合用PS-InSAR技術(shù)對(duì)高家灣滑坡區(qū)近幾年變形情況進(jìn)行監(jiān)測(cè)分析,而Envisat ASAR數(shù)據(jù)周期較長(zhǎng)且只有歷史時(shí)期的存檔數(shù)據(jù),用SABS-InSAR技術(shù)處理可以得到較理想的計(jì)算結(jié)果,更好地監(jiān)測(cè)歷史時(shí)期變形情況。數(shù)據(jù)基本特征如表1所示。

2.1 SENTINEL-1A數(shù)據(jù)

數(shù)據(jù)空間基線、多普勒中心和數(shù)據(jù)的時(shí)間分布見(jiàn)圖2,其中空間基線以0為中心,正態(tài)分布于兩側(cè),空間基線都在150m以?xún)?nèi),為圖像高精度配準(zhǔn)提供了基礎(chǔ);多普勒中心分布于0到-0.1之間;以2016年10月27日為參考主影像。

2.2 ENVISAT ASAR 數(shù)據(jù)

考慮到該區(qū)域和SAR數(shù)據(jù)的特點(diǎn),確保有效數(shù)據(jù)被用于構(gòu)建小基線集。最終選定空間基線閾值700米,時(shí)間基線545天,以2009年4月14日為參考主影像。所得數(shù)據(jù)集結(jié)構(gòu)較為合理,可以滿(mǎn)足SBAS技術(shù)要求。如圖3所示。

2.3 技術(shù)方法

如何提高相干性,有效消除各種誤差,成為高精度InSAR變形觀測(cè)的關(guān)鍵,對(duì)此,一是通過(guò)改善算法細(xì)節(jié)來(lái)實(shí)現(xiàn),二是通過(guò)構(gòu)建全新的框架和算法來(lái)提高相干性,獲取更詳細(xì)、準(zhǔn)確的變形信息。在構(gòu)建新的框架和算法方面,近年來(lái),主要是在傳統(tǒng)D-InSAR技術(shù)的基礎(chǔ)上進(jìn)一步發(fā)展了時(shí)序InSAR技術(shù),這類(lèi)技術(shù)通過(guò)長(zhǎng)時(shí)間序列的InSAR數(shù)據(jù)進(jìn)行分析,解決低相干性及去除大氣、DEM誤差的影響,極大推動(dòng)了InSAR技術(shù)的發(fā)展[10]。包括四個(gè)層面的算法:(1)永久散射體干涉測(cè)量及技術(shù)(Persistent Scatter Interferometry,PSI);(2)StaMPS永久散射體技術(shù)(Stanford Method for Persistent Scatterers);(3)SBAS(Small Baseline Subset)技術(shù);(4)PSI和SBAS相結(jié)合技術(shù)。

2.3.1 PS-InSAR技術(shù)

永久散射體干涉技術(shù)是意大利Ricca教授研究小組的Ferretti首先提出的[11]。PS-InSAR技術(shù)的基本原理是利用多景同一地區(qū)的SAR影像,通過(guò)統(tǒng)計(jì)分析時(shí)間序列上幅度和相位信息的穩(wěn)定性,探測(cè)不受時(shí)間、空間基線去相關(guān)影像的穩(wěn)定點(diǎn)目標(biāo)[12]。這些目標(biāo)點(diǎn)可能是人工建筑物、裸露的巖石、人工布設(shè)的角反射器等,由于它們?cè)跁r(shí)間序列SAR影像中幾乎不受斑點(diǎn)噪聲的影響,經(jīng)過(guò)很長(zhǎng)時(shí)間間隔仍然保持穩(wěn)定的散射特性,所以被稱(chēng)作PS(Persistent Scatter)。PS-InSAR可以精確估計(jì)并消除大氣效應(yīng)對(duì)相位的貢獻(xiàn),獲得米級(jí)精度的高程[13]和毫米級(jí)的形變量測(cè)[14]。永久散射體被定義為尺寸小于圖像像元尺寸,其散射系數(shù)強(qiáng)且長(zhǎng)期保持穩(wěn)定的地物目標(biāo)。其特點(diǎn)是在空間基線超過(guò)臨界基線的情況下,能保持不受空間失相干的影響,因此可利用所有可能獲得的圖像生成單主圖像堆棧(星星結(jié)構(gòu)),來(lái)獲取永久散射體0.5mm/a的地表變形信息[15]。

2.3.2 SBAS-InSAR技術(shù)

SBAS-InSAR是由Berardino[16]和Lanari[17]等研究人員提出的一種與PS-InSAR技術(shù)采用不同策略的時(shí)序InSAR分析方法,短基線技術(shù)采用小基線子集差分干涉測(cè)量技術(shù) [18];將所獲得的長(zhǎng)時(shí)間序列SAR數(shù)據(jù)組合成具有小空間基線干涉子集的(集合內(nèi)SAR圖像基線距小,集合間的SAR圖像基線距大)D-InSAR數(shù)據(jù)組進(jìn)行解纏計(jì)算,然后通過(guò)奇異值分解方法(SVD)將不同時(shí)間基線的空間小基線子集數(shù)據(jù)組合形成時(shí)間序列,進(jìn)而獲取高精度變形場(chǎng)。其針對(duì)的對(duì)象是面散射體(Distribution Scatter,DS),其精度較PS-InSAR低,但更適合于非城市區(qū)地表變形監(jiān)測(cè)。相比PS-InSAR技術(shù),SBAS方法只需要少量的SAR影像,而且由于生成干涉圖的影像間基線較短,能夠更好的避免空間失相干,正是由于研究區(qū)歷史SAR數(shù)據(jù)(EnviSAT ASAR)比較缺乏,因此對(duì)歷史變形過(guò)程的反演,SBAS方法對(duì)該區(qū)域有較好的適用性。

3 滑坡變形數(shù)據(jù)獲取

3.1 2003-2010年歷史形變

運(yùn)用SBAS-InSAR技術(shù)對(duì)27景ENVISAT ASAR數(shù)據(jù)進(jìn)行裁剪后,經(jīng)多次實(shí)驗(yàn)后選擇時(shí)間基線閾值為700天,空間基線不大于臨界基線的55%的小基線集網(wǎng)絡(luò)結(jié)構(gòu)較為合理。采用精密軌道數(shù)據(jù)進(jìn)行粗配準(zhǔn),90米分辨率的SRTM數(shù)據(jù)作為模擬地形相位的DEM,1:5的多視視數(shù),濾波窗口64*64,采用3D解纏算法,閾值0.35。經(jīng)過(guò)配準(zhǔn)、去平、濾波、軌道精煉、重去平、相位解纏和地理編碼等處理步驟,最終得到高家灣滑坡區(qū)歷史時(shí)期地表形變數(shù)據(jù)。

3.2 2014-2017年近期形變

使用PS-InSAR技術(shù)對(duì)57景sentinel-1A數(shù)據(jù)進(jìn)行處理,通過(guò)計(jì)算分析棄除大氣效應(yīng)的不同效果,選擇基于單主圖像的星型結(jié)構(gòu)圖像拓?fù)浣Y(jié)構(gòu)建立數(shù)據(jù)集最佳。參考DEM為90米分辨率的SRTM數(shù)據(jù),采用精確軌道參數(shù)進(jìn)行粗配準(zhǔn),選擇大于1500個(gè)的同名像素點(diǎn),建立多項(xiàng)式方程,進(jìn)行精配準(zhǔn),配準(zhǔn)精度達(dá)到1/8個(gè)像元。采用Modificated Goldstein方法,1*1的多視,15*15的濾波窗口。選擇星型拓?fù)浣Y(jié)構(gòu),稀疏點(diǎn)模式選擇Amp.Stab.Index模式,高程的閾值選擇-500米-500米區(qū)間,基于時(shí)間序列的濾波平滑的步長(zhǎng)選擇5,其中多項(xiàng)式的項(xiàng)數(shù)選擇1次多項(xiàng)式。

3.3 時(shí)間序列變形趨勢(shì)分析

采用PS-time統(tǒng)計(jì)工具計(jì)算分析高家灣滑坡區(qū)的地表形變數(shù)據(jù)PS點(diǎn)的變化趨勢(shì)。此方法是Matteo Berti 等人2013年[19]通過(guò)對(duì)在意大利北部亞平寧山脈收集的1000個(gè)時(shí)間序列樣本數(shù)據(jù)的統(tǒng)計(jì)分析,編寫(xiě)的一個(gè)應(yīng)用程序,在統(tǒng)計(jì)測(cè)試序列的基礎(chǔ)上提供永久散射體(PS點(diǎn))時(shí)間序列有條件的自動(dòng)分類(lèi)。程序有六個(gè)預(yù)定義的目標(biāo)趨勢(shì),即不相關(guān)、線性相關(guān)、二次相關(guān)、雙線性相關(guān)、勻速不連續(xù)和變速不連續(xù)。這些趨勢(shì)被認(rèn)為是PS點(diǎn)變形數(shù)據(jù)的“典型”模式。

4 監(jiān)測(cè)結(jié)果及分析

4.1 形變過(guò)程及趨勢(shì)分析

對(duì)于上述方法獲取的高家灣滑坡區(qū)的地表變形數(shù)據(jù),通過(guò)PS-time工具統(tǒng)計(jì)分析后,將研究區(qū)內(nèi)所有PS點(diǎn)的平均形變速率計(jì)算結(jié)果疊加顯示在高家灣滑坡區(qū)的DEM底圖上。結(jié)合野外勘測(cè)的滑坡方向、滑坡邊界和地裂縫等相關(guān)資料(由青海省國(guó)土資源廳提供),該區(qū)域內(nèi)的地表變形特征、變形分布規(guī)律和年均累積位移變化曲線如圖4示。

圖4顯示2004-2010年間,該滑坡區(qū)形變量0-3.67mm,比較輕微且分布分散,地表活動(dòng)相對(duì)比較穩(wěn)定。圖5顯示2014-2017年間滑坡區(qū)內(nèi)PS點(diǎn)速率變化值集中在-27.37-5.77mm,PS點(diǎn)各類(lèi)型變化趨勢(shì)總體呈遠(yuǎn)離雷達(dá)視線向的變形趨勢(shì)(下沉),與實(shí)地勘測(cè)結(jié)果比較吻合。滑坡邊界以外的PS點(diǎn)速率變化明顯比較小,多數(shù)在-2.8mm以下,說(shuō)明這些區(qū)域地表活動(dòng)相對(duì)穩(wěn)定,發(fā)生滑坡可能性較小。

圖6的統(tǒng)計(jì)分析處理結(jié)果顯示,2014年10月-2016年1月期間,累積位移走勢(shì)比較平緩,基本在-5mm以?xún)?nèi)小幅波動(dòng);2016年1月-2016年2月期間,累積位移走勢(shì)明顯加速,位移量從-5mm迅速下降到-12mm,此階段即地表變化比較劇烈的時(shí)期。推斷可能與2016年1月21日青海省海北州門(mén)源6.4級(jí)地震有關(guān)。2016年2月之后,累積位移又回歸平緩變化狀態(tài),總的位移量在-15-22mm間波動(dòng)。

4.2 主裂縫兩側(cè)不同部位變形特征

針對(duì)野外調(diào)查獲取的主裂縫情況,根據(jù)PS-InSAR結(jié)果,對(duì)主裂縫兩側(cè)分別提取PS點(diǎn)進(jìn)行分析,選點(diǎn)區(qū)域及PS點(diǎn)分布情況(圖7),得到主裂縫兩側(cè)突變點(diǎn)位移均值曲線(圖8)。

由位移曲線可以看出,裂縫兩側(cè)形變趨勢(shì)有較大差異,尤其是在2016年1月之后主裂縫兩側(cè)位移差持續(xù)增大,裂縫右側(cè)區(qū)域加速下沉,左側(cè)區(qū)域下沉速度相對(duì)較小。此時(shí)刻與野外調(diào)查結(jié)果吻合,與期間滑坡區(qū)累積位移走勢(shì)明顯加速,地表變化比較劇烈的特征吻合,與2016年1月21日青海省海北州門(mén)源6.4級(jí)地震的發(fā)生時(shí)間吻合。

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

時(shí)序InSAR技術(shù),能夠獲得緩慢的地表變形數(shù)據(jù),對(duì)滑坡監(jiān)測(cè)效果較好。通過(guò)多次對(duì)比實(shí)驗(yàn),不斷優(yōu)化干涉測(cè)量計(jì)算參數(shù)配置,能夠在高家灣滑坡區(qū)內(nèi)獲得較高精度和高相干性的永久散射點(diǎn)(PS點(diǎn)),作為該滑坡區(qū)地表形變監(jiān)測(cè)的數(shù)據(jù)基礎(chǔ)。通過(guò)PS-time統(tǒng)計(jì)分析后,可以較好歸納出滑坡區(qū)地表形變的時(shí)空特征和變化趨勢(shì),對(duì)滑坡監(jiān)測(cè)分析很有意義。

雖然利用InSAR技術(shù)獲取滑坡區(qū)地表的形變信息,結(jié)合野外調(diào)查資料進(jìn)行了比對(duì)。但是,衛(wèi)星獲取的形變方向?yàn)槔走_(dá)視線方向,衛(wèi)星的軌道固定,山體不同部分的坡度坡向不一致,導(dǎo)致獲取的形變結(jié)果會(huì)與坡體形變的真實(shí)情況存在一定的差異。限于沒(méi)有GPS野外觀測(cè)數(shù)據(jù),無(wú)法驗(yàn)證形變結(jié)果精度,得到的只是該滑坡區(qū)地表的相對(duì)形變。

參考文獻(xiàn):

[1]王治華.滑坡遙感[M].北京:科學(xué)出版社,2012.

[2]王立偉.基于D-InSAR數(shù)據(jù)分析的高山峽谷區(qū)域滑坡位移識(shí)別[D].北京科技大學(xué),2015.

[3]王志勇,張金枝.基于INSAR技術(shù)的滑坡災(zāi)害監(jiān)測(cè)[J].大地測(cè)量與地球動(dòng)力學(xué),2013,33(3):87-91.

[4]楊長(zhǎng)江,易祎,趙蓉.基于Sentinel-1數(shù)據(jù)的時(shí)序InSAR技術(shù)在滑坡監(jiān)測(cè)方面的應(yīng)用——以巴東地區(qū)為例[J].科技創(chuàng)新與生產(chǎn)力,2017(04):1674-9146.

[5]張?jiān)娗眩Y建軍,繆亞敏,等.基于SBAS技術(shù)的岷江流域潛在滑坡識(shí)別[J].山地學(xué)報(bào),2018,36(1):91-97.

[6]F Amelung, DL Galloway, JW Bell, HA Zebker, RJ Laczniak. Sensing the ups and downs of Las Vegas: InSAR reveals structural control of land subsidence and aquifer-system deformation. Geology, 1999,27(6):483-486.

[7]DA Schmidt, R Bürgmann. Time-dependent land uplift and subsidence in the Santa Clara valley, California, from a large interferometric synthetic aperture radar data set. Journal of Geophysical Research Solid Earth,2003,108(B9).

[8]ZW Li, XL Ding, GX Liu. Modeling atmospheric effects on InSAR with meteorological and continuous GPS observations: algorithms and some test results. Journal of Atmospheric and Solar-Terrestrial Ph..., 2004,66(11):907-917.

[9]HA Zebker, PA Rosen, S Hensley. Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps. Journal of Geophysical Research Solid Earth,1997,102(B4):7547-7563.

[10]許才軍,何平,溫?fù)P茂,等.InSAR技術(shù)及應(yīng)用研究進(jìn)展[J].測(cè)繪地理信息,2015,40(2):1-9.

[11]C Colesanti, A Ferretti, F Ferrucci, C Prati, F Rocca.Monitoring known seismic faults using the permanent scatterers (PS) technique. IEEE International Geoscience & Remote Sensing ..., 2000,5:2221-2223 vol.5.

[12]廖明生,王騰.時(shí)間序列INSAR技術(shù)與應(yīng)用[M].北京:科學(xué)出版社,2014.

[13]Perssin D.Validation of the sub-metric accurary of vertical positioning of PS's in C band. IEEE Transactions on Geoscience and Romete Sensing etters,2008,5(3):502-206.

[14]Ferreti A, Savio G, et al. Submillimeter Accuracy of InSAR Time Series: Experimental Validation: Experimental validation. IEEE Transactions on Geoscience and Remote Sensing,2007,45(5):1142-1153.

[15]A Ferretti, C Prati, F Rocca.Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience & Remote Sensing,2001,39(1):8-20.

[16]P Berardino, G Fornaro, R Lanari ,E Sansosti.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Transactions on Geoscience & Remote Sensing,2003,40 (11):2375-2383.

[17]R Lanari, O Mora, M Manunta, J Mallorquí, P Berardino, ...A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms. Geoscience & Remote Sensing IEEE Transactions on,2004,42(7):1377-1386.

[18]P Berardino, G Fornaro, R Lanari, E Sansosti.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Transactions on Geoscience & Remote Sensing,2003,40(11):2375-2383.

[19]M. Berti, A. Corsini S. Franceschini, and J. P. Iannacone .Automated classification of Persistent Scatterers Interferometry time series. Nature Hazards and Earth System Sciences, 2013,13(8):1945-1958.

猜你喜歡
變形
變形記
談詩(shī)的變形
柯西不等式的變形及應(yīng)用
“變形記”教你變形
不會(huì)變形的云
“我”的變形計(jì)
會(huì)變形的折紙
童話世界(2018年14期)2018-05-29 00:48:08
變形巧算
例談拼圖與整式變形
會(huì)變形的餅
主站蜘蛛池模板: 伊人久久大香线蕉影院| 欧美精品影院| 亚洲AV无码久久精品色欲| 欧美午夜一区| 国产18在线| 成人免费午夜视频| 精品国产免费第一区二区三区日韩| 久久久久无码国产精品不卡 | 婷婷综合缴情亚洲五月伊| 国产一区二区三区免费观看| 91网址在线播放| 久久免费精品琪琪| 国产农村精品一级毛片视频| 亚洲色图狠狠干| 国产视频a| 五月激激激综合网色播免费| 国产对白刺激真实精品91| 国产一区二区三区视频| 日韩欧美高清视频| 日本欧美成人免费| 日韩欧美高清视频| 四虎成人免费毛片| 高清乱码精品福利在线视频| 99激情网| 夜色爽爽影院18禁妓女影院| 欧美区一区| 久久精品国产精品青草app| 久久久久亚洲av成人网人人软件| 久久久亚洲色| 波多野结衣在线se| 在线色综合| 色有码无码视频| 色综合天天操| 国产 在线视频无码| 欧美国产精品不卡在线观看| 她的性爱视频| 成人av专区精品无码国产| 亚洲成A人V欧美综合| 制服丝袜在线视频香蕉| 国产精品尤物在线| 少妇被粗大的猛烈进出免费视频| 久久婷婷色综合老司机| 国产成人一区免费观看| 国产毛片高清一级国语| 亚洲二三区| 亚洲色图综合在线| 99国产精品一区二区| 亚洲一区二区无码视频| 夜夜操国产| 国产午夜福利片在线观看 | 亚洲码在线中文在线观看| 欧美成人国产| 老司国产精品视频91| 久久99精品国产麻豆宅宅| 国产成人亚洲综合a∨婷婷| 亚洲av日韩av制服丝袜| 男女男精品视频| 五月婷婷丁香色| 深爱婷婷激情网| 又黄又湿又爽的视频| 国产精品极品美女自在线看免费一区二区| 日本高清成本人视频一区| 动漫精品啪啪一区二区三区| 香蕉视频在线观看www| 波多野结衣在线一区二区| AV无码无在线观看免费| 国产精品无码制服丝袜| 成年人国产视频| 四虎永久在线视频| 亚洲最大福利视频网| 狠狠v日韩v欧美v| 欧美一区日韩一区中文字幕页| 毛片久久网站小视频| 找国产毛片看| 夜精品a一区二区三区| 国产一区亚洲一区| 五月婷婷综合网| 久久亚洲综合伊人| 在线网站18禁| 午夜啪啪网| 国产二级毛片| 久久久久青草大香线综合精品 |