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

基于不同建模數據的陜西省志丹縣滑坡易發性分區

2021-01-16 02:49:09王麗英王紅梅郭盈盈紀丁愈
人民長江 2021年12期
關鍵詞:模型研究

王麗英 王紅梅 郭盈盈 紀丁愈

摘要:滑坡易發性分區是預測滑坡的有效方法,但目前的建模數據較為單一,而建模數據會對分區的結果造成影響,因此,需要對建模數據源進行深入挖掘。以陜西省延安市志丹縣為研究區,選擇了10種滑坡誘發因子,基于289個滑坡樣本,分別計算了各誘發因子的頻率密度(FR)和盒維數,構建了2組差異化的建模數據;隨后利用核函數邏輯回歸模型(KLR)制作了研究區滑坡易發性分區圖;最后采用平均絕對誤差(MAE)和ROC曲線下的面積(AUC)對分區的結果和模型進行了評價和對比。結果表明:利用盒維數作為建模數據可有效提升分區結果的精度和模型的分類以及泛化能力,值得在研究區推廣。

關 鍵 詞:滑坡; 易發性分區; ROC曲線; 誘發因子; 頻率密度; 盒維數; 志丹縣; 陜西省

中圖法分類號: P694

文獻標志碼: A

DOI:10.16232/j.cnki.1001-4179.2021.12.015

0 引 言

黃土滑坡是陜西省北部的黃土丘陵溝壑區范圍內最具威脅的自然災害之一,并且由于黃土本身的結構特性以及黃土高原特殊的自然環境,滑坡發生概率較高[1]。志丹縣位于陜西省延安市西北部,地貌類型屬于典型的黃土丘陵溝壑區,全縣范圍內地質災害頻發,根據志丹縣2019年地質災害防治方案顯示,2019年內列入全縣地質災害群測群防數據庫的隱患點共39處,其中滑坡22處,在每年在汛期時,都會發育大量的滑坡災害,嚴重威脅當地居民的生產生活,制約著當地經濟的發展。并且當地人口集中,加之資源和土地的需求量大,而地質災害發育的范圍仍在擴張,因此開展滑坡災害的易發性分區工作十分必要。

對于滑坡易發性分區的研究,國內外學者不僅充分發揮了GIS技術在可視化分析、圖形編輯、數據管理、空間分析、DEM模型等功能方面的優勢,還先后提出了信息量模型[2]、專家打分模型[3]、邏輯回歸模型[4]、判別分析模型[5]、人工神經網絡模型[6]、支持向量機模型[7]、物元模型[8]等,基于不同的尺度開展了不同地區的滑坡易發性分區工作。雖然這些研究取得了豐碩的成果,但是在滑坡易發性分區建模的樣本準備階段,建模數據是影響易發性分區結果的關鍵因素[9],而現有的研究所使用的建模數據較為單一,還未對不同類型的建模數據進行深入的挖掘和對比。

鑒于此,本文以志丹縣作為研究區,分別基于兩組差異化的建模數據集,構建核函數邏輯回歸模型(Kernel logistic regression model,KLR),繪制研究區的滑坡易發性分區圖。在評估易發性分區精度的基礎上,對比評價建模數據對分區結果的影響。所得到的結果可以為當地地質災害防治工作提供建議,同時也可以為黃土丘陵溝壑區的滑坡易發性分區研究提供一種新的方法。

1 研究區概況和數據來源

1.1 研究區概況

志丹縣位于陜西省延安市西北部的黃土丘陵溝壑地區,地理坐標位于東經108°11′56″~109°3′48″,北緯 36°21′23″~37°11′47″之間(見圖1)。全縣東西長約70 km,南北寬約92.56 km,總面積約3 781 km2。志丹縣的氣候類型屬于典型的溫帶大陸性季風氣候,年平均氣溫為8.1 ℃,最低氣溫為-28 ℃,最高氣溫為37.4 ℃。由于志丹縣地處黃土高原,年平均降雨量約為472.2 mm,降雨的月份在一年之中分布不均勻,6~8月為降雨相對集中的月份,降雨量約為292.9 mm,降雨通常以暴雨形式出現,占年降雨量的56%。降雨區域分布不均勻,總體趨勢為從西北部地區到東南部地區降雨量逐漸增多,愈偏北降雨量愈少。流經志丹縣境內的河流主要有洛河、周河、杏子河,縣內以洛河、周河、杏子河3條河流為主干,大小支流、沖溝極為發育,溝壑密度為1.3 km/ km2,河流比降為1.6‰。

志丹縣地貌類型屬于以梁峁為主體的黃土梁峁丘陵溝壑區地貌,在長期的侵蝕作用下形成了幾個顯著的特點:地表破碎、梁峁密布、溝壑縱橫、河谷深切、基巖裸露。全縣范圍內主要發育中生代-新生代地層,具體包括三疊系、侏羅系、新近系和第四系地層,其中第四系黃土覆蓋于全縣境內,其余時代的底層大多零星的出露于沖溝或者河谷兩側。根據地層年代將研究區的巖性類型劃分為4組(見表1)。志丹縣的地質構造類型屬于華北陸臺鄂爾多斯地臺中的陜北盆地,構造類型是以延安地區為中心的陜北單斜翹曲構造,整體形狀呈東高西低的趨勢。據《陜西省區域地質志》記載[10],志丹縣境內無4級以上地震活動發生,所以本文研究的滑坡均為降雨型滑坡。

1.2 數據來源

本文所使用的原始數據包括:① 3.24 m×3.24 m分辨率的GF-2全色多光譜遙感影像(軌道號921309,日期2018-07-29);②30 m×30 m分辨率的數字高程模型(DEM);③ 1∶100 000比例尺的研究區土地利用圖;④1∶200 000比例尺的研究區地質圖;⑤ 基于野外實際調查的289個滑坡樣本。

1.3 樣本預處理

滑坡在遙感影像上所呈現出的平面形態是不規則的,由于受到研究區圖幅比例尺的影響,當研究范圍大于1 000 km2時,平面面積最大的滑坡在圖中也無法顯現出其全部的輪廓。并且在建模運算時,柵格圖形分辨率的高低也會對模型運算的速度造成影響[11]。因此,對于滑坡易發性分區而言,有必要在建模前期對樣本進行預處理。研究區內發育的滑坡總面積為30.24 km2,約占研究區總面積的0.8%。因此,為了提高運算效率,本文采用質心法(Centroid method)將研究區內289個不規則的滑坡圖斑轉換為289個滑坡點。

由于后續的滑坡易發性分區涉及到因子篩選和機器學習建模,而建立和評估模型的重點就是要將數據分別劃分為訓練樣本集(Training dataset)和測試樣本集(Validation dataset),所以需要創建與正樣本(滑坡點)同等數量的負樣本(非滑坡點),并按照7/3的比例對滑坡點和非滑坡點進行隨機分割,最終得到由404個樣本(202個滑坡點和202個非滑坡點)構成的訓練樣本集和由174個樣本(87個滑坡點和87個非滑坡點)構成的測試樣本集(見圖1)。

1.4 誘發因子提取與分級

基于研究區的地質環境背景資料,選取了高程、平面曲率、剖面曲率、坡度、坡向、道路緩沖區、水系緩沖區、年平均降雨量、NDVI、地層巖性作為研究區滑坡誘發因子。隨后利用ArcGIS軟件,基于所選用的數據源對誘發因子進行提取,并將因子圖層全部重采樣為30 m×30 m分辨率(見圖2)。最后對提取出的誘發因子進行分級。

2 研究方法

2.1 FR的計算方法

利用頻率密度(Frequency ratio,FR)量化滑坡誘發因子是一種較為常用的方法。

通常以Sim/S表示誘發因子各分級內的滑坡樣本數(Sim)占總樣本數(S)的百分比,以Zim/Z表示誘發因子各分級的柵格數(Zim)占總柵格數(Z)的百分比,利用公式(1) 可計算得出各誘發因子分級對應的FR。

FR=Sim/SZim/Z(1)

隨后,基于計算得到的FR值,分別對每一種滑坡誘發因子進行重分類,構建出基于FR的建模數據集1。

2.2 分維數的測量方法

分維數是定量反映功能、形態、時間、信息、空間等特征在局部與局部、局部與整體之間具有統計意義的參數,也可以定量反映滑坡的空間分布特征[12]。為了測量研究區滑坡在空間內分布的分形維數,本文采用目前應用最廣泛的分維數測量方法,即盒維數法(Box-counting Method)。盒維數法是既可以在二維、三維空間中進行分維數估算的經典方法,同時也適用于點集[13]。

盒維數的測量原理為:首先利用邊長為γ的正方形對分割研究區,統計包含滑坡點的網格數N(γ),不斷改變γ對研究區重新劃分,統計有滑坡點分布的網格數量,從而獲得點對序列(γ,N(γ))。按照1 m的間隔改變γ 50次,使得點對序列滿足或者近似滿足公式(2),則可計算出對應的盒維數d。

Nγ=γ-d(2)

同樣的,基于測算出的盒維數,分別對每一種滑坡誘發因子進行重分類,構建出基于盒維數的建模數據集2。

2.3 誘發因子優選

在滑坡易發性分區建模的過程中,并不是所有的滑坡誘發因子都對模型有積極的貢獻[14]。因此,為了提高模型運算的精度和效率,有必要對誘發因子進行篩選。本文通過計算各誘發因子對應的信息增益率(Information gain ratio,IGR)來完成因子篩選。

熵可以反映滑坡誘發因子在評價體系中的混亂程度,而信息增益率是通過構建C4.5決策樹,定量測算各誘發因子的熵增速率,從而判斷誘發因子對滑坡易發性分區模型的貢獻程度。IGR的輸出范圍為[0,1],當誘發因子的IGR為0時,表明該因子對模型無貢獻,需要排除,且不參與建模。

2.4 KLR模型簡介

KLR模型的原理是在傳統邏輯回歸模型的基礎上,利用核函數(Kernel Function)將研究區內的滑坡點轉換至一個n維空間內,使得線性不可分的問題得以在高維空間中被解決。KLR模型的計算方法如下:

LogitP=w·φx+b(3)

這里災損土地發生的可能性同樣用P表示,w和b是模型的兩個參數,φ為核函數,則公式(3) 的對數形式可以寫為

P=11+exp(w·φx+b)(4)

本文所使用的核函數為RBF核函數,其表達式如下:

Kxi,xj=exp-δxi-xj2,δ>0(5)

這里δ控制著RBF核函數的敏感度,在計算時需要手動輸入。

2.5 模型對比和結果評估的方法

2.5.1 統計學指標

滑坡易發性分區結果的精度會對后續滑坡的防治工作造成影響,因此有必要對分區的結果進行評估。本文利用平均絕對誤差(Mean Absolute Error,MAE)來對易發性分區的精度進行檢測[15]。MAE表示所有獨立樣本的觀測值與目標值之間的偏差,可以定量反映實際預測誤差的大小,MAE越小,說明結果的精度越高。

MAE=1nnz=1Az-F(z)(6)

式中:n為滑坡樣本點的數量,z表示滑坡樣本點,A(z)和F(z)分別表示觀測值和目標值。

2.5.2 ROC曲線

本文利用受試者接收特征曲線(Receiver Operating Characteristic Curve,ROC)對滑坡易發性分區模型進行比較。ROC曲線在滑坡易發性分區研究中應用十分廣泛,同時也取得了令人滿意的效果[16]。ROC曲線的橫軸和縱軸分別為1-特異度(1-Specificity)和敏感度(Sensitivity),通常利用曲線下面積(AUC)來判定模型的優劣[17]。AUC的計算方法如下:

Sensitivty=TPTP+FN(7)

Specificity=TNTN+FP(8)

AUC=TP+TNTP+TN+FP+FN(9)

式中:TP和TN分別表示被正確分類的滑坡點和非滑坡點的數量,FN和FP分別表示被錯誤分類的滑坡點和非滑坡點的數量。一般來講,AUC越接近1,說明模型的分類能力越強。

3 結果與分析

3.1 因子篩選結果

基于數據集1和數據集2中的訓練數據,分別計算數據集中各滑坡誘發因子對應的IGR。從表2可以看出,在數據集1和數據集2中,道路緩沖區因子對模型的貢獻度都最高,IGR值分別為0.688 3和0.635 8;數據集1中得平面曲率因子對模型的貢獻度最低(IGR為0.115 4),數據集2中的剖面曲率因子對模型的貢獻度最低(IGR為0.107 7)。此外,研究區內所有滑坡誘發因子的IGR值均大于0,且對應的標準差均處于合理的范圍內,表明所有的誘發因子都對模型有積極的貢獻。因此,保留數據集1和數據集2中的全部滑坡誘發因子。

3.2 基于差異化建模數據的滑坡易發性分區

基于表3所列出的FR以及盒維數,構建了數據集1和數據集2。分別利用數據集1和數據集2中的訓練數據訓練KLR模型。為了對模型中的參數進行優調,本文采用10折交叉驗證(10-Cross Validation)的方法分別獲得了數據集1對應的參數(δ=0.025 2)和數據集2對應的參數(δ=0.027 1)。

分別利用數據集1和數據集2構建KLR模型,模型輸出的后驗概率即為滑坡易發性指數(Landslide susceptibility index,LSI),輸出區間為[0,1]。LSI越趨近于1,說明滑坡發生的概率越高,反之亦然。隨后基于自然間斷點法將兩組LSI分割為5個區間,分別代表極低易發區、低易發區、中易發區、高易發區和極高易發區。最后利用ArcGIS軟件對劃分后的易發區進行可視化,結果如圖3所示。

3.3 分區結果精度評價

本文分別利用數據集1和數據集2中的訓練數據和測試數據,計算分區結果對應的MAE來完成精度評價。從訓練數據檢測的結果可以看出(見圖4):基于數據集1所構建的KLR模型的MAE=0.24,基于數據集2所構建的KLR模型的MAE=0.18。從測試數據檢測的結果可以看出(見圖5):基于數據集1所構建的KLR模型的MAE=0.31,基于數據集2所構建的KLR模型的MAE=0.11。結果表明,利用數據集2作為建模數據生成的滑坡易發性分區結果的精度高于數據集1。

3.4 模型對比

本文利用ROC曲線下的面積(AUC)對比模型的分類能力和泛化性。分別基于數據集1和數據集2中的訓練數據繪制ROC曲線,見圖6(a),結果顯示數據集2對應的AUC值最高為0.880 6,并且95% CI和標準差均最低,表明利用數據集2構建的KLR模型對滑坡的分類能力最強。

分別基于數據集1和數據集2中的測試數據繪制ROC曲線,見圖6(b),結果顯示數據集2對應的AUC值同樣最高為0.920 3,且95% CI和標準差均最低,表明利用數據集2構建的KLR模型的泛化性最優,值得在研究區推廣。

4 討 論

在滑坡易發性分區建模前,排除無貢獻程度的因子是非常必要的[18]。眾多研究結果表明:如果利用對模型無貢獻的誘發因子參與建模,不僅會降低運算效率,增加過擬合的風險,同時也會對分區結果的精度造成影響[19-20]。因此本文通過計算信息增益率來定量評估滑坡誘發因子對模型的貢獻程度,排除無貢獻的誘發因子,以保證分區結果的精度。

KLR模型是基于LR模型衍生出的一種高性能機器學習算法,具有結構簡單等優點。KLR模型不僅在滑坡災害的研究中被廣泛使用,同時也在其他的研究領域被使用[21]。在構建KLR模型的過程中,參數的選擇直接決定著分類的精度[22]。因此本文利用了10折交叉驗證的方法進行調參,目的就是為了將參數對模型的影響降至最低。另外,選擇不同的核函數也會產生不同的結果,鑒于此,在今后的研究中,會嘗試利用不同的核函數來完成滑坡易發性分區,并對結果進行評估。

利用FR作為模型的建模數據,在滑坡研究中十分常見,但是利用盒維數作為KLR模型的建模數據的研究還少之又少。從張庭瑜等[23-24]的研究可以看出,在縣域尺度下,盒維數不僅可以有效地反映出滑坡在空間中的分布狀態,同時也可以揭示滑坡誘發因子與滑坡之間的內在聯系。另外從結果可以明顯看出,利用盒維數作為模型的建模數據,可以取得更高精度的分區結果,同時,模型的分類能力和泛化性也優于傳統建模數據。因此,認為利用盒維數作為滑坡易發性分區的建模數據是一種高效且新穎的方法。

5 結 論

本文分別利用頻率密度和盒維數量化滑坡誘發因子,構建了2種差異化的建模數據集(數據集1和數據集2)。以核函數邏輯回歸模型(KLR)和289個滑坡樣本為基礎,分別利用兩種數據集構建了KLR模型,制作了陜西省延安市志丹縣滑坡易發性分區圖。利用平均絕對誤差(MAE)和ROC曲線下的面積(AUC)對分區的結果和模型進行了評價。結果表明:當利用盒維數作為建模數據時,滑坡易發性分區結果的平均絕對誤差最小(MAE=0.11),精度最高,且模型對滑坡的分類能力最強,泛化性最優。同時所得到的結果可以為當地的滑坡防治工作提供數據基礎,也可以為今后研究區內的滑坡研究提供參考。

參考文獻:

[1] 徐張建,林在貫,張茂省.中國黃土與黃土滑坡[J].巖石力學與工程學報,2007,17(7):1297-1312.

[2] 張向營,張春山,孟華君,等.基于GIS和信息量模型的京張高鐵滑坡易發性評價[J].地質力學學報,2018,34(1):96-105.

[3] 張向營,張春山,孟華君,等.基于Random Forest和AHP的貴德縣北部山區滑坡危險性評價[J].水文地質工程地質,2018,45(4):142-149.

[4] TSANGARATOS P,ILIA I.Comparison of a logistic regression and Nave Bayes classifier in landslide susceptibility assessments:the influence of models complexity and training dataset size[J].Catena,2016,145(2):164-179.

[5] 李秀珍,王成華,宋剛.基于Fisher判別分析法的潛在滑坡判識模型及其應用[J].中國地質災害與防治學報,2009,20(4):23-26.

[6] CHEN W,POURGHASEMI H R,ZHAO Z.A GIS-based comparative study of Dempster-Shafer,logistic regression and artificial neural network models for landslide susceptibility mapping[J].Geocarto International,2017,32(4):367-385.

[7] KUMAR D,THAKUR M,DUBEY C S,et al.Landslide susceptibility mapping & prediction using Support Vector Machine for Mandakini River Basin,Garhwal Himalaya,India[J].Geomorphology,2017,295(15):115-125.

[8] 吳益平,殷坤龍,姜瑋.浙江省永嘉縣滑坡災害風險預警研究[J].自然災害學報,2009,18(2):124-130.

[9] JUNG-HYUN L,SAMEEN M I,PRADHAN B,et al.Modeling landslide susceptibility in data-scarce environments using optimized data mining and statistical methods[J].Geomorphology,2018,303(9):284-298.

[10] 陜西省地質礦產局.陜西省區域地質志[M].北京:地質出版社,1989.

[11] HYUN-JOO O,SARO L,SOO-MIN H.Landslide susceptibility assessment using frequency ratio technique with Iterative random sampling[J].Journal of Sensors,2017,2017(4):1-21.

[12] YANG Z Y,POURGHASEMI H R,LEE Y H,et al.Fractal analysis of rainfall-induced landslide and debris flow spread distribution in the Chenyulan Creek Basin,Taiwan[J].Journal of Earth Science,2016,27(1):151-159.

[13] LIU L,LI S,LI X,et al.An integrated approach for landslide susceptibility mapping by considering spatial correlation and fractal distribution of clustered landslide data[J].Landslides,2019,16(4):1-14.

[14] HONG H,PRADHAN B,SAMEEN M I,et al.Improving the accuracy of landslide susceptibility model using a novel region-partitioning approach[J].Landslides,2018,15(4):753-772.

[15] PHAM B T,PRAKASH I,SINGH S K,et al.Landslide susceptibility modeling using Reduced Error Pruning Trees and different ensemble techniques:hybrid machine learning approaches[J].Catena,2019,175(5):203-218.

[16] 韓玲,張庭瑜,張恒.基于IOE和SVM模型的府谷鎮滑坡易發性分區[J].水土保持研究,2019,26(3):373-378.

[17] 張庭瑜,韓玲,張恒,等.混合分類模型在滑坡易發性分區中的適用性研究:以延安市寶塔區為例[J].干旱區資源與環境,2020,22(1):192-201.

[18] YOUSSEF,AHMED M,AL-KATHERY,et al.Landslide susceptibility mapping at Al-Hasher area,Jizan(Saudi Arabia)using GIS-based frequency ratio and index of entropy models[J].Geoscience Journal,2015,19(1):113-134.

[19] KAVZOGLU T,KUTLUG SAHIN E,COLKESEN I.An assessment of multivariate and bivariate approaches in landslide susceptibility mapping:a case study of Duzkoy district[J].Natural Hazards,2015,76(1):471-496.

[20] YOUSSEF A M,POURGHASEMI H R,POURTAGHI Z S,et al.Landslide susceptibility mapping using random forest,boosted regression tree,classification and regression tree,and general linear models and comparison of their performance at Wadi Tayyah Basin,Asir Region,Saudi Arabia[J].Landslides,2016,13(5):839-856.

[21] CHEN W,XIE X,PENG J,et al.GIS-based landslide susceptibility modelling:a comparative assessment of kernel logistic regression,Nave-Bayes tree,and alternating decision tree models[J].Geomatics Natural Hazards & Risk,2017,8(2):950-973.

[22] WANG L J,SAWADA K,MORIGUCHI S.Landslide-susceptibility analysis using light detection and ranging-derived digital elevation models and logistic regression models:a case study in Mizunami City,Japan[J].Journal of Applied Remote Sensing,2013,7(1):35-61.

[23] ZHANG T Y,HAN L,HAN J,et al.Assessment of landslide susceptibility using integrated ensemble fractal dimension with Kernel Logistic Regression Model[J].Entropy,2019,21(2):218-234.

[24] ZHANG T Y,HAN L,ZHANG H,et al.GIS-based landslide susceptibility mapping using hybrid integration approaches of fractal dimension with index of entropy and support vector machine[J].Journal of Mountain Science,2019,16(6):1275-1288.

(編輯:劉 媛)

Landslide susceptibility mapping based on differentiated modeling data in Zhidan County,Shaanxi Province

WANG Liying1,WANG Hongmei2,GUO Yingying1,JI Dingyu3

(1.School of Intelligent Construction,Chongqing Vocational College of Civil Engineering,Chongqing 400072,China; 2.Urban Construction College,Chongqing Energy Vocational College,Chongqing 402260,China; 3.State Key Laboratory of Hydraulics and Mountain River Development and Protection,Sichuan University,Chengdu 610065,China)

Abstract:

Landslide susceptibility mapping is an effective method to predict landslides,but the current modeling data is relatively single,moreover the modeling data may affect the results of the mapping.Therefore,it is necessary to dig deeper into the modeling data sources.In view of this,we selected 10 landslide triggering factors in Zhidan County,Yanan City,Shaanxi Province.Based on 289 landslide samples,the frequency ratio and box counting dimension of each triggering factor were calculated,and two sets of differentiated modeling data were constructed.Then,the landslide susceptibility maps of the study area were made by using the kernel function logistic regression model.Finally,the average absolute error (MAE) and the area under the ROC curve (AUC) were used to evaluate and compare the results of landslide susceptibility mapping.The results show that using the box dimension as the modeling input data can effectively improve the accuracy of the partition results and the classification and generalization capabilities of the model,which is worthy of promotion in the study area.

Key words:

landslide;susceptibility mapping;ROC curve;landslide triggering factors;frequency ratio;box counting dimension;Zhidan County;Shaanxi Province

猜你喜歡
模型研究
一半模型
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
新版C-NCAP側面碰撞假人損傷研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产制服丝袜91在线| 亚洲a免费| 欧美视频在线第一页| 久久黄色视频影| 天堂岛国av无码免费无禁网站| 91 九色视频丝袜| 久久综合伊人 六十路| 精品人妻系列无码专区久久| 人妖无码第一页| 亚洲欧美日韩另类在线一| 日韩精品成人在线| 成人毛片在线播放| 色婷婷电影网| 九九免费观看全部免费视频| 一区二区三区高清视频国产女人| 国产区人妖精品人妖精品视频| 国产精品3p视频| 亚洲男人的天堂网| 精品国产黑色丝袜高跟鞋| 国产无码精品在线| 久久久精品无码一区二区三区| 欧美第二区| 久久中文字幕2021精品| 免费无码又爽又黄又刺激网站| 国产成人永久免费视频| 国产亚洲欧美在线专区| 无码专区国产精品第一页| 2020极品精品国产| 2020久久国产综合精品swag| 怡红院美国分院一区二区| 97久久精品人人| 99re在线免费视频| 久久精品中文字幕免费| 亚洲,国产,日韩,综合一区 | 欧美激情视频一区| 无遮挡国产高潮视频免费观看| 美女视频黄又黄又免费高清| 伦精品一区二区三区视频| 国产不卡在线看| 日本高清有码人妻| 国产一级毛片高清完整视频版| 久久精品电影| 亚洲Av激情网五月天| 三级毛片在线播放| 九九九国产| Aⅴ无码专区在线观看| 国内精品伊人久久久久7777人| 亚洲免费成人网| 青青草国产精品久久久久| 91精品啪在线观看国产91九色| 在线亚洲小视频| 国产国语一级毛片| 在线五月婷婷| 青青草原国产| 日本一区高清| 亚洲天堂精品视频| 久久鸭综合久久国产| 亚洲日本一本dvd高清| 日韩人妻无码制服丝袜视频| 欧美成人精品欧美一级乱黄| 国产精品黑色丝袜的老师| 国内嫩模私拍精品视频| 亚洲精品国产综合99久久夜夜嗨| 在线精品欧美日韩| 成人综合在线观看| 国内精品久久久久久久久久影视 | 激情乱人伦| 欧美国产日韩另类| 香蕉99国内自产自拍视频| 一区二区偷拍美女撒尿视频| 中国丰满人妻无码束缚啪啪| 高清国产va日韩亚洲免费午夜电影| 欧洲一区二区三区无码| 色偷偷一区| 欧美精品啪啪一区二区三区| 在线观看网站国产| 嫩草国产在线| 免费观看精品视频999| 一级毛片免费的| 草草线在成年免费视频2| 国产99久久亚洲综合精品西瓜tv| 久久综合丝袜长腿丝袜|