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

基于證據權和卡方自動交互檢測決策樹的滑坡易發性預測

2022-06-21 01:57:12黃發明石雨歐陽慰平洪安宇曾子強徐富剛
土木與環境工程學報 2022年5期
關鍵詞:關聯環境模型

黃發明,石雨,歐陽慰平,洪安宇,曾子強,徐富剛

(南昌大學 建筑工程學院,南昌 330031)

如何有效開展滑坡易發性預測制圖是現階段全世界范圍內區域滑坡研究的重點和難點。通過將GIS與數據驅動模型相結合,以圖像和數字的方式可構建出更高效準確的易發性預測模型。該易發性制圖的思路對滑坡高發地區的防災減災規劃具有重要意義[1]。

滑坡易發性可定義為特定地點在環境因子非線性耦合作用下發生滑坡的空間概率。基于地理相似性規律,即“地理環境越相似,地理特征越相近”可知,通過已經發生滑坡的環境因子來建立預測模型,則潛在滑坡的空間位置有可能被預測[2]。很明顯,從滑坡樣本點中確定滑坡易發性與其環境因子的關系式是易發性預測的關鍵所在,因此,選擇用以獲取輸入變量的滑坡-環境因子關聯分析法非常重要。隨著遙感和GIS等基礎數據源獲取技術的進步,易發性建模的空間數據源及其質量有了較大提升[3]。一般而言,具體研究區內的滑坡環境因子類型可通過相關文獻綜述和研究區的自然地理和地質條件確定。筆者重點關注滑坡易發性建模過程中滑坡與其環境因子的非線性關聯分析這一不確定性因素,并進一步研究其對滑坡易發性建模的影響。

啟發式模型、數理統計模型和機器學習模型是易發性預測過程中常用的3種類型[4]。啟發式模型[5]和數理統計模型被大量使用,主要有確定性因子(Certainty Factors,CF)[6]、層次分析法[7]和多元線性回歸[5]等;機器學習相關模型包括邏輯回歸(Logistic Regression,LR)[8]、C5.0決策樹[9]、人工神經網絡[10-11]、隨機森林(Random Forest,RF)[12]、支持向量機(Support Vector Machines,SVM)[13]、卡方自動交互檢測(Chi-squared Automatic Interaction Detector,CHAID)決策樹[14]和貝葉斯網絡[15]等。對于哪種類型的模型最適合易發性預測,現階段還沒有一致的意見,但優秀的機器學習模型能夠提高滑坡易發性預測精度,對滑坡易發性區間劃分有著顯著影響,并可能進一步改變滑坡易發性級別的劃分。筆者擬用CHAID決策樹這一被廣泛應用的典型機器學習方法構建滑坡易發性模型并探索建模不確定性特征。

在將建模預測出的滑坡易發性指數(Landslide Susceptibility Index,LSIs)與各類環境因子開展聯系時,需開展滑坡與其基礎環境因子(不考慮誘發因子)之間的非線性關聯分析,其關聯值可直接作為易發性模型的輸入變量[16]。目前,常用的關聯分析法包括確定系數[17]、頻率比(Frequency Ratio,FR)[18]、熵指數(Index of Entropy,IOE)[16]和證據權重(Weight of Evidence,WOE)[19]等。不同關聯分析法的內部計算思路具有較大的差異性,導致各方法下的易發性建模存在不確定性[20-21]。關聯分析法太粗糙會導致部分信息丟失,降低模型預測精度;優秀的關聯分析法能獲取較準確的環境因子影響滑坡發育的信息,進一步提高滑坡環境因子分析及其建模的可靠性。可見,探討不同關聯分析法對易發性預測建模的影響規律具有重要意義。

學者們采用不同關聯分析法和模型開展易發性預測建模,例如:Zhang等[22]應用IOE模型、LR-IOE和SVM-IOE模型獲得了中國陜西省府谷縣滑坡易發性圖,結果表明,LR-IOE模型的準確率最高,其次是IOE模型和SVM-IOE模型。李文彬等[23]深入探討滑坡與其環境因子間的非線性聯接以及不同數據驅動模型對滑坡易發性預測建模不確定性的影響規律,結果表明,RF模型預測性能最優,WOE-RF模型預測的滑坡易發性不確定性較低。張鐘遠等[24]基于地理信息系統平臺構建了云南省鎮康縣滑坡易發性預測指標體系,結果顯示,頻率比耦合LR模型具有更高的成功率和預測率。但大多數情況下,現有研究使用特定的關聯分析法開展易發性預測建模,而較少提供可信的依據和合理的解釋,并且較少深入探討這種不確定性因素對易發性預測建模的影響。通過探討關聯分析法耦合模型下的滑坡易發性結果的不確定性,更能深入理解易發性預測的可靠性和可行性,可降低關聯分析法不確定性因素帶來的影響。

筆者采用FR和WOE兩種非線性關聯分析法的計算數據值與原始環境因子數據(以下簡稱“原始因子數據”)作為CHAID決策樹模型的輸入變量,以陜西省延長縣為例,開展滑坡易發性預測建模的不確定性分析,包括精度評價、LSIs分布規律和平均秩等。

1 滑坡易發性建模分析

FR和WOE兩種關聯法耦合CHAID決策樹模型時的易發性預測建模流程(圖1)如下:

1)獲取研究區滑坡編錄及相關環境因子數據源以便構建易發性建模的空間數據集;

2)將FR、WOE和原始因子數據作為CHAID決策樹的輸入變量,形成3種耦合模型;

3)分別對3種耦合模型開展易發性預測建模,然后在GIS中繪制滑坡易發性圖并劃分易發性等級;

4)通過ROC精度、均值、標準差和平均秩等對易發性預測結果進行不確定分析;

5)通過對比分析找到最佳關聯分析法,為易發性建模提供指導。

圖1 滑坡易發性預測建模流程圖Fig.1 Flowchart of landslide susceptibility prediction

1.1 滑坡與環境因子的關聯分析法

1.1.1 頻率比 頻率比(Frequency Ratio,FR)反映了滑坡在各環境因子類別的分布狀況,闡述環境因子各屬性區間對滑坡的相對影響度,并且能夠很好地解釋滑坡與各因子之間的內在聯系[25]。FR>1代表在對應的環境因子條件下利于滑坡事件的發生;FR<1表明該環境因子區間的屬性與滑坡的發展關系較弱。利用環境因子的FR值作為各模型的輸入變量之一,其計算公式如式(1)。

(1)

式中:Nj為環境因子某區間中出現的滑坡柵格數;N是全區已知滑坡所分布柵格的總數;Sj是環境因子的單元數;S是全區柵格總數。

1.1.2 證據權 證據權(Weight of Evidence,WOE)法在貝葉斯準則基礎上綜合各類證據層來實現定量計算某事件的發生概率。WOE法通過將滑坡編錄和各類環境因子層進行空間關聯,從而得到滑坡處各環境因子的詳細分布特征權重因子W+和W-,其在每個環境因子分級中的計算如式(2)、式(3)所示。

(2)

(3)

1.2 卡方自動交互檢測決策樹

CHAID決策樹以卡方統計量為基礎實現最優決策樹構建,也就是通過自變量和因變量間的解釋性來實現因變量的自動判別。CHAID決策樹具有強大的非線性擬合預測性能,能容忍樣本數據缺失及樣本量不足等缺陷。CHAID模型設定樹生長的層數、分裂及聚合閾值等停止標準來構建準確高效的預測或分類模型,同時,為防止過擬合現象而用隨機分成的訓練樣本構建模型;最后再利用隨機分成的測試樣本對CHAID進行逐步檢驗,以修正模型參數。

1.3 不確定性分析方法

1.3.1 ROC 曲線精度分析 采用受試者工作特征(Receiver Operating Characteristic,ROC)曲線下面積(Area Under ROC,AUC)值作為一種量化指標來整體評估建模性能。ROC曲線對測試集中各樣例進行排序并依序選擇各截斷點,再逐個把樣例作為正例來進行計算,依據當前分類器的“真陽率”和“假陽率”進行ROC曲線的繪制,相關評價指標如表1所示。AUC值等于隨機挑選的正樣本的排名高于隨機挑選的負樣本的概率,AUC值越大,則易發性模型預測性能越好[4]。

表1 ROC曲線的相關指標

1.3.2 易發性指數統計規律分析 均值(Mean)是集中趨勢的測量,計算如式(4)所示(式中:Xn為第n個柵格單元的滑坡易發性指數值),其量化了研究區LSIs分布的整體偏向趨勢,反映了LSIs分布的平均水平。標準差(Standard Deviation)是對圍繞平均值的離差的測量,計算如式(5)所示(式中:μ為滑坡易發性指數均值;Xi為第i個柵格單元的滑坡易發性指數值),量化了LSIs分布的離散程度,標準差越小,說明LSIs越接近平均值,反之,則說明其與平均值的差異越大。采用均值和標準差從整體上分析LSIs的分布特征,揭示不同關聯分析法和模型耦合模型下的預測性能,為滑坡易發性研究提供理論指導[23]。

(4)

(5)

1.3.3 易發性指數的差異顯著性 采用顯著性差異水平進一步分析各耦合模型下易發性建模的不確定性。具體采用Kendall協同系數檢驗法,對任意兩組不同耦合模型下預測出的LSIs進行差異顯著性檢驗。若Kendall秩相關系數W小于1及檢驗結果的顯著性小于0.05,說明這兩組耦合模型下LSIs的差異是顯著的,拒絕原假設。本文通過成對因子顯著性檢驗發現,W值為0.139,小于1,且P值均小于0.05,可見,各耦合模型下的LSIs間差異顯著[27]。

2 延長縣簡介及環境因子分析

2.1 延長縣簡介及滑坡編錄

延長縣位于陜西東部,面積約2 368.7 km2,地勢從西北向東南方向傾斜。縣境內屬黃土高原丘陵溝壑區(河谷階地、黃土溝谷區、黃土溝間區和巖質丘陵區),出露三疊系中上統內陸湖相碎屑沉積巖和第四系風積、沖洪積和堆積黃土等地層,新近系砂礫巖在研究區出露較少(圖2)。另外,縣境內地質構造活動強度低,屬于暖溫帶干旱大陸性季風氣候,年均降雨量約564 mm且集中在7、8、9月份。

圖2 延長縣滑坡編錄圖Fig.2 Landslide inventory map of Yanchang

根據已有的滑坡野外調查資料和數據庫可知,延長縣共發生滑坡82處,主要類型為小型淺層覆蓋滑坡,主要運動方式為牽引式(59%)和推移式滑動(41%);縣境內的小型滑坡45處(占比54.8%),中型滑坡36處(占比43.9%),大型滑坡只有1處。延長縣滑坡分布位置如圖2所示,滑坡主要分布在縣域西部及周邊地區,東部和中部較少;大部分發生滑坡的位置地勢較高,距離河流水系也較近。延長縣滑坡的發生與地層巖性和工程活動密切相關。

2.2 環境因子分析

2.2.1 環境因子介紹 根據延長境內滑坡的特征及相關參考文獻的介紹,利用遙感影像和GIS軟件系統從數據源中提取14類滑坡環境因子,包括地形、水文、地表覆被和基礎地質等[28-29]。其中,高程、NDVI、NDBI和MNDWI等12個因子為連續型數據,而距河流距離和地層巖性2個因子為離散型數據(表2)。對于連續型環境因子,先通過小間隔對該因子進行等分,再依據FR和WOE值將數值相近的區間合并成一個類別[30]。對于離散型數據類型的環境因子,采用固有的自然分組來進行分級:距河流距離因子按照距河流距離100、300、400、500、800、900、1 000 m和大于1 000 m進行分類;地層巖性因子為三疊系砂巖夾砂質泥巖和油頁巖(T2t)、三疊系厚層砂巖夾泥巖(T3h)、三疊系細砂層粉砂巖夾與泥巖互層(T3y)、三疊系厚層狀長石石英砂巖(T2w)和第四系更新統風積和洪積黃土(Qp1-3)[31]。另外,在使用原始因子數據作為CHAID決策樹模型的輸入變量時,將距河流的距離和地層巖性兩種離散型數據類型的環境因子進行了“啞變量”處理。

2.2.2 地形地貌因子 高程、坡度、坡向、剖面曲率、平面曲率、地形起伏度、地形粗糙度、地形切割深度和地形濕度指數等環境因子均從DEM中提取(圖3)[23,32]。以地形起伏度為例,分析其8個等級區間內的FR和WOE值(表2),發現滑坡發生概率與研究區的地形起伏度大小成正比。在20~4區間內發生滑坡的概率最大,為78.34%;其中,FR值均大于1,WOE值均為正值,35~40區域內FR和WOE值最大,分別為2.843和1.148。FR和WOE值都顯示出地形起伏度大小與滑坡發生有著較強的正向相關性,可見關聯分析法在表達滑坡與地形起伏度的非線性關聯性時具有較為一致的趨勢和計算效果。

表2 環境因子的關聯分析值

續表2

續表2

圖3 延長縣滑坡環境因子Fig.3 Landslide environmental factors in Yanchang

2.2.3 水文環境因子 由于河流對邊坡的浸潤和侵蝕作用,越靠近河流的邊坡土壤含水量可能越高,導致斜坡體失穩的可能性更高[33-34]。利用距河流距離和MNDWI來表征水文環境對滑坡發育的影響。以距河流的距離因子為例(表2),當距河流距離小于400 m時,滑坡發育的可能性更高(達74.41%),其中,FR值均大于1,WOE值均為正值;在100~300 m區域內,FR和WOE值最大,分別為1.873和0.992。

2.2.4 地表覆被因子 NDBI和NDVI分別反映了研究區域內的建筑分布和自然植被對滑坡地質災害發育的影響[35]。從表2可知,當NDVI在0.121~0.424范圍內時,其與滑坡有較強的關系,該區間包括了研究區內近年來所有的已發生的滑坡;其中,在0.121~0.182范圍內,FR值大于1且WOE值為正數。NDBI能較好地反映研究區域內建筑的分布情況,當NBVI在0.730~0.949范圍內時幾乎囊括了近年來研究區內所有的滑坡,間接反映了人類工程建設對滑坡發育的影響。

2.2.5 基礎地質因子 巖土類型表征滑坡體的物質基礎[36-37],分析表2可知,T3h和T3y巖性區域面積僅占延長縣面積的10.6%,而區域內滑坡發生的概率高達23.2%,且FR值均大于1、WOE值均為正值,說明T3h和T3y巖性區域內滑坡發生的頻率較高;在Qp1-3巖性條件下,滑坡發生概率高達76.8%;在T2t巖性區域內,無滑坡分布;T2w巖性區域在研究區內占比比較小,結果不具有研究意義。

3 延長縣滑坡易發性預測建模

3.1 數據準備

30 m分辨率的柵格被廣泛用作滑坡易發性的制圖單元,基于30 m分辨率,整個延長縣被劃分為2 622 482個柵格,已發生的82處滑坡被劃分為3 403個滑坡柵格[38]。通過FR和WOE兩種關聯法對14個環境因子各屬性區間進行重新賦值,作為CHAID決策樹開展易發性建模的輸入變量;同時,也以原始因子數據作為輸入變量開展單獨CHAID決策樹的滑坡易發性建模。通過SPSS modeler 18.0軟件把3 403個滑坡柵格單元賦值為1,同時隨機挑選與滑坡單元相同數量的非滑坡單元,并將其易發性賦值為0,作為模型輸出變量;然后按7∶3隨機劃分滑坡和非滑坡柵格單元(6 806個)及其相關屬性值,得到模型訓練集和測試集。最后將整個研究區柵格單元的FR和WOE關聯分析值以及原始因子數據代入訓練好的模型中,預測延長縣LSIs,并將其按照自然間斷點法[39]劃分為5個易發性級別。

3.2 延長滑坡易發性預測結果

在SPSS modeler軟件中進行CHAID決策樹建模。以WOE樣本數據為例,首先需從外部源中讀取源節點,將6 806個滑坡-非滑坡樣本數據導入SPSS modeler軟件中;接著對字段屬性、測量級別及各字段在建模中的角色進行選擇或修改;再經由分區選擇將樣本數據分為訓練集(70%)和測試集(30%);然后在CHAID建模節點字段選項卡中使用預定義角色,應用boosting算法創建一個整體,由其生成模型序列以增強模型預測的準確度;選擇CHAID樹生長算法并定制樹的最大深度值為5、父節點的最小記錄數為75、子節點的最小記錄數為15,以此來限制決策樹的增長;CHAID決策樹的其他參數使用SPSS modeler中的默認值;最后將整體環境因子的WOE帶入訓練好的CHAID決策樹模型中,實現延長縣滑坡LSIs的準確預測。FR-CHAID和單獨CHAID決策樹模型的建模步驟和參數設置與WOE-CHAID決策樹模型基本一致。

3.3 滑坡易發性制圖表達

分兩步開展滑坡易發性制圖,首先將3種耦合模型預測出的LSIs導入GIS軟件中,然后依據自然間斷點法將延長縣滑坡易發性劃分為極高、高、中等、低和極低5類等級區間[33]。WOE-CHAID、FR-CHAID和單獨CHAID決策樹模型下的滑坡易發性結果如圖4所示。延長縣大部分地區屬于低和極低易發區,滑坡高和極高易發區主要位于坡度和高程中等且距離河流較近的山地丘陵地區。但3種耦合模型下得到的滑坡易發性級別存在顯著差異,圖4中延長縣內已發生的82處滑坡幾乎都落在WOE-CHAID和FR-CHAID決策樹模型預測的極高與高易發性等級區域內,而單獨CHAID決策樹模型預測的極高與高易發性等級區域與82處滑坡位置存在些許偏差。

圖5 CHAID決策樹模型的滑坡易發性建模ROC曲線Fig.5 ROC curve of landslide susceptibility modeling of CHAID decision tree

4 滑坡易發性預測不確定性分析

4.1 ROC精度評價

采用測試集AUC值作為具體指標量化不同耦合模型的預測性能,AUC值越大,表明耦合模型預測性能越優。WOE-CHAID、FR-CHAID和單獨CHAID決策樹模型的滑坡易發性結果ROC曲線如圖5所示。從圖5中可知,3種耦合模型下的結果均較好且相對穩定,表現出良好的滑坡易發性性能。AUC精度從大到小依次為:AUC(WOE-CHAID)>AUC(FR-CHAID)>AUC(單獨的CHAID),說明FR和WOE兩種關聯分析法在CHAID決策樹模型中具有比原始因子數據更穩定的易發性預測性能。WOE耦合CHAID決策樹模型的易發性預測效果最好且預測效率最高,AUC精度較FR提高了2.1%,較原始因子數據提高了3.1%。

4.2 滑坡易發性指數分布規律

采用均值和標準差分別反映LSIs分布的平均水平和離散程度,并以此分析耦合模型下的易發性預測不確定性。WOE-CHAID、FR-CHAID和單獨CHAID決策樹模型預測的LSIs分布不確定性規律較為一致,在極低和低易發區分布較集中而在高和極高易發區分布逐漸減少。LSIs平均值從小到大排名為:單獨的CHAID (0.364)

圖6 CHAID決策樹模型的LSIs分布Fig.6 LSIs distribution of CHAID decision tree

4.3 耦合模型預測易發性指數的差異性分析

采用顯著性差異水平來進一步分析各耦合模型下易發性建模的不確定性,通過該試驗計算各耦合模型下預測的LSIs的平均秩,以便對易發性模型性能排序。平均秩越小則模型性能越好,最終模型比較結果為:WOE-CHAID決策樹模型預測LSIs的平均秩(值為1.85)最小,其次是FR-CHAID(值為2.06) 和單獨的CHAID決策樹(值為2.09)模型。顯著性差異水平和平均秩顯示出各耦合模型的易發性建模存在不確定性,如何規避這些不確定性是獲得可靠的易發性模型的重要研究內容。

4.4 滑坡環境因子重要性分析

滑坡環境因子的重要性反映了已發生的滑坡事件受該環境因子影響程度的大小[40]。由于原始因子數據和不同的關聯分析值在易發性預測建模中有著不同的表現,基于CHAID決策樹模型中自帶的分類器屬性來評估在原始因子數據、FR和WOE等輸入變量下各個環境因子的重要性。另外,易發性建模中共使用14個環境因子(原始因子數據含“啞變量”類型,共23個環境因子),排名10名之后的環境因子重要性均小于0.04,因此僅展示重要性排名前10的環境因子。從圖7可知,坡度、地形起伏度、距河流的距離(原始因子數據中為100~300 m和500~800 m的兩個“啞變量”因子)、地形切割深度和地形粗糙度等5個環境因子在單獨CHAID、FR-CHAID和WOE-CHAID決策樹易發性預測中有著較大的貢獻,占據重要性排名均在前5位,重要性均大于0.08。其次,平面曲率和地形濕度指數在所有決策樹模型中也發揮著相對重要的作用,重要性均大于0.04。

圖7 滑坡環境因子重要性Fig.7 The importance of environmental factors of

4.5 各關聯分析法的性能分析

關聯分析法通過定量統計可直觀表現各環境因子不同屬性區間對滑坡易發性空間的影響性。Li等[27]、Saha等[41]對上述部分關聯分析法反映滑坡與其環境因子空間關聯的性能進行了對比分析,所得結果與筆者研究基本一致。由上述分析可知,環境因子與滑坡間的空間信息的關聯性表達越充分,則LSIs的區分度越大,進一步的易發性預測效果就越佳。在FR和WOE關聯分析法的環境因子分級中,WOE更能反映環境因子內部影響滑坡發育的空間信息的差異,具有更優的預測精度(AUC=86.3%);FR相較于WOE法更加簡潔高效,在保證易發性精度的同時能有效避免太復雜的統計分析;基于原始因子數據進行的單獨CHAID決策樹模型易發性預測精度略小于FR-CHAID和WOE-CHAID決策樹模型。此外,單獨的CHAID、FR-CHAID和WOE-CHAID決策樹模型預測的LSIs平均值逐漸減小而標準差逐漸增大,且平均秩也逐漸減小。可見關聯分析法的易發性預測建模效果較好,WOE優于FR,而原始因子數據的易發性建模效果較差。

由文獻[27,42]可知,滑坡與環境因子(不考慮誘發因子)之間的非線性關聯分析法種類繁多。筆者僅使用FR和WOE兩種關聯分析法耦合CHAID決策樹模型進行滑坡易發性的不確定性對比分析而并未考慮其他關聯分析法,在下一步研究中可以考慮使用概率法、信息量、確定性系數和熵指數等其他關聯分析法,耦合多種不同類型的模型開展更加全面的易發性預測不確定性分析。

5 結論

1)WOE-CHAID決策樹模型易發性預測的AUC精度最高,且均值和平均秩較小,標準差較大;FR-CHAID決策樹的AUC精度略低于WOE-CHAID,可見WOE具有更優秀的非線性關聯性能。

2)將原始因子直接用作輸入變量的單獨CHAID決策樹模型的易發性預測精度整體略低于關聯分析法的耦合模型。為了提高滑坡易發性建模效率,可直接使用單獨CHAID決策樹模型,但要體現滑坡與其環境因子的空間關聯性或分析環境因子各子區間對滑坡發育的影響規律,則使用關聯分析法和CHAID決策樹模型耦合建模的優勢顯著。

3)總體來說,WOE-CHAID決策樹模型的易發性預測結果可靠性最高,預測出的LSIs與實際的滑坡概率分布特征更加相符。

猜你喜歡
關聯環境模型
一半模型
“苦”的關聯
當代陜西(2021年17期)2021-11-06 03:21:36
長期鍛煉創造體內抑癌環境
一種用于自主學習的虛擬仿真環境
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
孕期遠離容易致畸的環境
環境
奇趣搭配
智趣
讀者(2017年5期)2017-02-15 18:04:18
主站蜘蛛池模板: 日本一区二区三区精品AⅤ| 黄片在线永久| 国产人成午夜免费看| 国内老司机精品视频在线播出| 久久国产精品电影| 欧美国产精品不卡在线观看| 国产福利在线观看精品| 毛片免费在线视频| 国产91蝌蚪窝| 精品偷拍一区二区| 欧美精品1区2区| 亚洲最新地址| 中文字幕资源站| 国产免费福利网站| 美女一区二区在线观看| 国产女人在线观看| AV在线麻免费观看网站| 久久99蜜桃精品久久久久小说| 久久久久九九精品影院| 视频一区亚洲| 经典三级久久| 国产av无码日韩av无码网站| 无码一区二区三区视频在线播放| 热这里只有精品国产热门精品| 国产内射一区亚洲| 久久综合亚洲鲁鲁九月天| 国产呦精品一区二区三区下载 | 色综合五月婷婷| 白丝美女办公室高潮喷水视频 | 就去色综合| 亚洲第一综合天堂另类专| 天堂av综合网| 国产无码网站在线观看| 青青久视频| 日本91视频| 免费啪啪网址| 国产免费观看av大片的网站| 国产欧美日韩在线在线不卡视频| 成人av手机在线观看| 久久永久免费人妻精品| 欧美午夜在线视频| 美女国产在线| 欧美在线视频不卡第一页| 美女国内精品自产拍在线播放 | 国产在线拍偷自揄拍精品| 成人国内精品久久久久影院| 久草视频一区| 在线视频亚洲色图| 尤物午夜福利视频| 国产三级国产精品国产普男人| 玩两个丰满老熟女久久网| 欧美国产日产一区二区| 美女无遮挡拍拍拍免费视频| 米奇精品一区二区三区| 国产嫖妓91东北老熟女久久一| 精品国产欧美精品v| 欧美人人干| 精品在线免费播放| 无码区日韩专区免费系列| 动漫精品啪啪一区二区三区| 天天色综网| 午夜福利在线观看入口| 国产99视频在线| 国产黑丝一区| 伊人网址在线| 99热这里只有精品在线播放| 久久青草精品一区二区三区| 国产日韩欧美在线播放| 无码专区国产精品一区| 亚洲午夜片| a级高清毛片| 四虎国产精品永久一区| 国产在线专区| 国产欧美另类| 欧美激情网址| 中文字幕 日韩 欧美| 精品国产美女福到在线不卡f| 天天操天天噜| 国产精品久久久久久久久| 国产精女同一区二区三区久| 波多野结衣中文字幕一区| 国产一国产一有一级毛片视频|