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

基于隨機森林模型的安徽省土壤屬性空間分布預測①

2019-07-26 02:57:58盧宏亮趙明松劉斌寅陸龍妹
土壤 2019年3期
關鍵詞:模型

盧宏亮,趙明松,2,劉斌寅,張 平,陸龍妹

基于隨機森林模型的安徽省土壤屬性空間分布預測①

盧宏亮1,趙明松1,2*,劉斌寅1,張 平1,陸龍妹1

(1 安徽理工大學測繪學院,安徽淮南 232001;2 土壤與農業(yè)可持續(xù)發(fā)展國家重點實驗室(中國科學院南京土壤研究所),南京 210008)

為探討隨機森林(random forest,RF)模型對土壤屬性空間預測的精度,本文以安徽省為例,收集140個土壤樣本,利用GIS和RS技術,獲取相關的地形因子、遙感植被指數(shù)及氣候數(shù)據(jù),利用RF模型分析土壤有機碳(SOC)含量、土壤容重和土壤黏粒含量與地形因子、遙感植被指數(shù)及氣候數(shù)據(jù)之間的關系,并進行空間分布預測。研究結果表明:①RF建模預測中,當節(jié)點分裂次數(shù)(mtry)值為1,決策樹數(shù)量(ntree)值分別為100、1 000和100時,獲得的SOC含量、土壤容重和土壤黏粒含量RF模型最優(yōu);②高程、歸一化植被指數(shù)(NDVI)、地貌、多尺度山谷平坦指數(shù)(MrVBF)和土壤類型是SOC含量的重要預測因子;地貌、年均降水量(MAP)、MrVBF、高程和土壤類型是土壤容重的重要預測因子;高程、MAP、MrVBF和平面曲率是土壤黏粒含量的重要預測因子;③RF模型可以較好地進行土壤屬性空間預測,多源環(huán)境變量組合可以分別解釋SOC含量、土壤容重和土壤黏粒含量的26%、23% 和22%;同時RF模型對于土壤類型和地貌等類型變量的處理具有一定優(yōu)勢。研究表明,在大尺度研究區(qū)域內,利用RF模型進行土壤屬性空間預測有一定的意義。

土壤屬性預測;隨機森林模型;環(huán)境變量;安徽省

土壤有機碳(soil organic carbon,SOC)是陸地生態(tài)系統(tǒng)平衡的主要因子[1],研究SOC含量的空間分布及其影響因素是陸地生態(tài)系統(tǒng)碳循環(huán)的基礎。土壤容重(bulk density,BD)是土壤的基本物理性質之一,對土壤的透氣性、入滲性能、持水能力、溶質遷移特征以及土壤的抗侵蝕能力有重要影響。土壤黏粒是土壤中最活躍的礦物成分[2],研究不同階段的土壤黏粒含量可以得出可靠的土壤相對年齡。研究上述土壤屬性的空間變異及其分布特征和環(huán)境因子的關系,對于了解生態(tài)系統(tǒng)、制定農業(yè)政策、進行土壤管理和監(jiān)測由于土地利用導致的環(huán)境變化有重要意義。

基于機器學習方法預測土壤屬性的空間分布逐步成為近年來的研究熱點。文獻研究表明,已有研究中使用的主要機器學習技術有分類和回歸樹(classification and regression tree,CART)[3]、k最近鄰(k-nearest neighbors,KNN)算法[4]、樹模型[5]和隨機森林(random forest,RF)模型[6]等:RF模型與大多數(shù)統(tǒng)計建模方法相比具有一些優(yōu)勢,它具有對多元共線性不敏感和不易出現(xiàn)過擬合問題等特點[6],且在噪聲和數(shù)據(jù)簡化處理方面最準確和穩(wěn)定[7]。國外研究中,Dharumarajan等[8]利用RF模型對印度南部半干旱熱帶地區(qū)的SOC、土壤pH等屬性進行了預測,結果證明RF模型可以提高土壤屬性空間預測的精度。Chagas等[9]比較了RF模型和多元線性回歸方法在半干旱地區(qū)土壤質地的空間預測制圖的效果,結果表明RF模型可以避免過擬合且預測精度更高。國內的研究中,郭彭濤等[10]基于多源環(huán)境變量和RF模型預測了橡膠園土壤全氮含量的空間分布,結果證明RF模型相較于逐步回歸、廣義加性混合模型和分類回歸樹等模型具有更高的預測性。姜賽平等[11]比較了普通克里格、回歸克里格、RF等模型在海南島土壤有機質(SOM)的預測研究中的精度,結果表明RF和回歸克里格模型能夠更好地描述SOM的局部變異信息。RF模型最合適進行土壤屬性空間的預測。研究RF模型在土壤屬性空間預測中的應用對數(shù)字土壤制圖具有一定的意義。

安徽省氣候差異明顯、地貌類型眾多、土地利用存在明顯的區(qū)域差異,這些條件的組合導致了多種環(huán)境因子共同影響土壤屬性的空間分布及變異。本研究以安徽省為例,利用GIS和RS技術提取土壤景觀環(huán)境因子,通過收集土壤野外采樣數(shù)據(jù),以SOC含量、土壤容重和土壤黏粒含量為預測目標,運用方差分析和相關性分析研究環(huán)境變量與預測目標的關系,通過RF建模選擇最優(yōu)環(huán)境變量組合和模型參數(shù)建立土壤屬性的預測模型并進行空間分布預測,同時探討省域尺度上3種土壤屬性的主要影響因素。研究結果有望為安徽省生態(tài)系統(tǒng)研究、土壤質量管理和農業(yè)生產提供數(shù)據(jù)基礎。

1 材料與方法

1.1 區(qū)域概況

安徽省(114°54' ~ 19°37'E,29°41' ~ 34°38'N)地處我國東部,跨長江、淮河中下游,東臨以上海為中心的長江三角洲經濟區(qū),西接中原腹地。安徽省總面積13.96萬km2,其中農田69%,低山丘陵14%,湖泊17%。全省處于亞熱帶向溫帶過渡帶。年均氣溫14 ~ 16 ℃,年均降雨量800 ~ 1 800 mm。除安徽西南和南部丘陵地區(qū)外,海拔一般不超過100 m。安徽省從北到南分為淮河中游平原、江淮丘陵崗地、沿江平原、皖西大別山區(qū)、皖南丘陵地區(qū)等5個地理區(qū)域。主要的土壤類型有:水稻土、潮土、砂姜黑土、黃棕壤、黃褐土、棕壤、紅壤、黃壤、紫色土、石質土、粗骨土、石灰土、山地草甸土等。

1.2 數(shù)據(jù)來源

本研究數(shù)據(jù)包括野外土壤調查數(shù)據(jù)、氣候數(shù)據(jù)、地形數(shù)據(jù)和遙感植被指數(shù)。土壤調查數(shù)據(jù)來源于《中國土系志?安徽卷》[12],該調查數(shù)據(jù)按照隨機性、均勻性和代表性的原則在安徽省全省范圍采集典型土壤剖面,數(shù)據(jù)集包含采樣點位置、景觀條件和土壤理化性質,采樣時間為2010—2011年。本研究選擇140個樣點的表層土壤屬性為預測的目標變量,包括SOC含量、土壤容重及土壤黏粒含量。氣候數(shù)據(jù)主要包括:年均溫(MAT)、年均降水量(MAP)。數(shù)據(jù)來自中國農業(yè)科學院農業(yè)資源與農業(yè)區(qū)劃研究所中國生態(tài)環(huán)境背景層面建造項目完成的柵格數(shù)據(jù)(1 km分辨率),為1980—1999年的逐月平均值計算生成。在ArcGIS支持下,從上述環(huán)境變量的柵格數(shù)據(jù)中提取各樣點的相應環(huán)境屬性。地形數(shù)據(jù)來源于地理數(shù)據(jù)空間云(http://www.gscloud.cn)的SRTM數(shù)字高程模型(DEM),空間分辨率為90 m。利用ArcGIS 10.2提取坡向、坡度、高程、平面曲率和剖面曲率;利用SAGA GIS 6.3.0提取多尺度山谷平坦指數(shù)(MrVBF)、多尺度脊頂平坦指數(shù)(MrRTF)、地形濕度指數(shù)(TWI)及地形位置指數(shù)(TPI),其中坡向數(shù)據(jù)由DEM數(shù)據(jù)產品中的SRTMTPI 坡位產品直接提取。歸一化植被指數(shù)(NDVI)和增強植被指數(shù)(EVI)來源于MODIS陸地產品16 d合成植被指數(shù)(MOD13Q1),空間分辨率為250 m,時間為2010年8月。所有環(huán)境變量及土壤屬性空間預測結果,分辨率統(tǒng)一為250 m。

1.3 隨機森林模型

隨機森林(random forest,RF)模型具有提高預測精度、減少過擬合、對缺失數(shù)據(jù)和多元共線性不敏感,且具有簡單處理大量的定量和定性數(shù)據(jù)能力的優(yōu)點[13]。對于土壤類型和地貌等類型變量,多數(shù)回歸模型處理方式比較復雜,一些研究甚至找不到適合的定性指標進行定量化描述[14],在R軟件中編程建立的RF模型只需將定性變量轉為因子(factor) 直接用于模型即可。

本研究使用R語言中的Random Forest 4.6軟件包進行建模。140個樣點按8:2分為建模集和驗證集。RF模型采用boostrap的方法對于樣本進行放回抽樣。沒有被抽取的記錄會自動生成一個對照集,所以不需要進行交叉驗證[6]。RF模型中的兩個可調參數(shù)決策樹數(shù)量(ntree)和節(jié)點分裂次數(shù)(mtry)決定了模型的配置。

1.4 數(shù)據(jù)處理與分析

利用SPSS 22 for windows進行方差分析和相關性分析,研究環(huán)境變量對于土壤屬性的影響[15]。對于定量環(huán)境變量,在R軟件中,使用scale() 函數(shù)進行歸一化處理后,利用SPSS進行相關性分析,將MAP、MAT、地貌和土壤類型進行方差分析,其中MAP分為 <800 mm、800 ~ 900 mm、900 ~ 1 000 mm、1 000 ~ 1 100 mm、>1 100 mm 5個降雨帶,MAT分8 ~ 10 ℃、10 ~ 12 ℃、12 ~ 14 ℃、14 ~ 16 ℃、>16 ℃ 5個溫度帶,地貌數(shù)據(jù)根據(jù)安徽省地貌類型矢量圖在ArcGIS中提取至采樣點,根據(jù)5個地貌區(qū)域劃分,在SPSS中以1 ~ 5重新編碼后進行方差分析,土壤類型數(shù)據(jù)根據(jù)已有土壤數(shù)據(jù)在SPSS中以1 ~ 13重新編碼后進行方差分析。

1.5 精度評價

模型精度評價選用均方根誤差(RMSE)、平均絕對誤差(MAE)[16]以及決定系數(shù)(2)[17]3個指標,其中MAE和RMSE越小表明預測精度越高,建模集2用于評價建模的擬合精度,驗證集2用于評價預測精度及模型泛化能力。計算方法如下:

2 結果與討論

2.1 土壤屬性統(tǒng)計特征

表1為安徽省土壤屬性統(tǒng)計結果。SOC含量介于1.33 ~ 33.53 g/kg,平均含量為14.57 g/kg;中等變異,變異系數(shù)為52.06%。土壤容重范圍為0.59 ~ 1.56 g/cm3,中等變異,變異系數(shù)為11.68%。土壤黏粒含量范圍為42.73 ~ 552.76 g/kg,中等變異,變異系數(shù)為47.43%。

表1 安徽省土壤屬性基本統(tǒng)計特征

注:SOC:土壤有機碳;BD:土壤容重;Clay:土壤黏粒,下表同。

2.2 土壤屬性的影響因素分析

相關性分析結果(表2)表明,坡向、高程、MrRTF和MAP與SOC含量顯著相關(<0.05);土壤容重與NDVI、坡向、高程、TPI、MrVBF、MrRTF、MAP和MAT顯著相關(<0.05);土壤黏粒含量則和NDVI、EVI、坡度、MrRTF、MrVBF和MAP都具有顯著相關性(<0.05)。

表2 土壤屬性與環(huán)境因子的相關性分析

注:*表示相關性達到<0.05顯著水平,**表示相關性達到<0.01顯著水平(雙尾)。

方差分析結果(表3)表明[18],不同的MAT、MAP和土壤類型的SOC含量和容重均存在顯著差異(<0.05),其他各因子對兩者的變異性均有顯著影響,地貌對于SOC含量的變異性沒有顯著影響(= 0.18);對于土壤黏粒含量,因子影響均顯著。對于SOC含量,土壤類型的值最大,說明土壤類型對SOC含量和容重的變異性影響最大;影響土壤黏粒含量變異性的最重要因素為MAP。

表3 安徽省各因子影響土壤屬性的方差分析

2.3 環(huán)境變量的篩選及重要性分析

利用RF模型進行變量重要性排序,對重要性較低的環(huán)境變量進行排除后重復建模,選取最優(yōu)環(huán)境變量組合用于預測。最終確定高程和NDVI等8個環(huán)境變量作為自變量預測SOC含量;地貌、MAP和土壤類型等9個環(huán)境因子作為土壤容重的預測因子;高程和MAP等8個環(huán)境因子用于土壤黏粒含量的預測。相關性分析和方差分析結果表明,NDVI及地貌與SOC含量的相關性并不顯著,RF模型的重要性分析卻表明NDVI和地貌是影響SOC含量重要的環(huán)境因素,這是由于RF模型對多元共線性不敏感。在進行土壤屬性預測變量篩選時,應該結合土壤學專業(yè)知識選取。

預測因子重要性排序表明(圖1),RF模型以增長均方誤差(increased in mean squared error,IncMSE)為變量重要性衡量指標,該值越大則變量重要性最高。對于SOC含量,高程、NDVI和地貌等為主要影響因子,高程和NDVI影響最大。有研究表明[8],NDVI與SOC含量呈極顯著正相關關系,即該指數(shù)越大,SOC含量越高,所以NDVI在SOC預測模型中為主要影響因子。影響土壤容重的環(huán)境因子中,地貌為最主要的影響因子。容重主要受土壤質地、結構的影響,不同地貌的土壤質地和結構區(qū)別顯著,所以在影響因子重要性排序中地貌為首要因素。對于土壤黏粒含量,高程、MAP、MrVBF和平面曲率是主要的預測因子。

2.4 隨機森林模型參數(shù)設定

本文通過逐次試驗,確定RF模型中mtry和ntree參數(shù)的最優(yōu)值[19]。固定mtry值(分別設為1、2和3),逐次調整ntree值(分別設為100、500和1 000),進行3組9次試驗。為避免過擬合問題,本文通過比較建模集和驗證集的2,選擇兩者最為接近的結果作為最優(yōu)預測模型。結果表明(表4),當mtry值為1,ntree值為100時,SOC預測模型的建模集和驗證集2最為接近,表明此時的模型穩(wěn)定性最好;當mtry值為1,ntree值分別為1 000和100時,容重和土壤黏粒含量的預測模型最為穩(wěn)健。

圖1 RF模型中土壤屬性預測因子的重要性排序

表4 RF模型中節(jié)點分裂次數(shù)(mtry)和決策樹數(shù)量(ntree)的篩選

2.5 預測精度分析及空間分布

RF模型的性能通過計算RMSE、MAE、2等參數(shù)來進行評估,經過參數(shù)調優(yōu)后采用最穩(wěn)定的RF模型作為最終預測模型。結果(表5、圖2)表明:①驗證集中SOC含量、容重和黏粒含量的決定系數(shù)分別為0.27、0.22和0.21。建模集中的決定系數(shù)與驗證集相近,說明RF模型有效避免了過擬合的問題,這與前人的理論一致[6];②SOC含量預測效果最好,土壤容重次之,對土壤黏粒的預測效果最差;③對于SOC含量,建模集的2和驗證集的2均高于0.25且整體水平相近,說明模型擬合度和泛化能力均較高,且模型較穩(wěn)定;對于容重和土壤黏粒含量,建模集的2和驗證集的2基本相同,說明模型穩(wěn)定性極高,但是預測精度較低。④由MAE和RMSE可以看出,模型整體預測精度較高,說明在大尺度區(qū)域上,RF模型對于土壤屬性仍然有不錯的預測效果。

利用RF模型分別對安徽省SOC含量、容重和土壤黏粒含量進行預測得到三者的空間分布圖(圖3),其中圖A ~ C為實測值圖,D ~ F為預測值圖。由圖3可知,安徽省SOC含量分布大致為由北向南逐漸增加,這基本符合以往的研究[20-21],其中淮河中游平原地區(qū)SOC含量最低,沿江平原東部SOC含量最高。淮河中游平原地區(qū)土壤容重值最高,其他區(qū)域大致由北向南逐漸降低。安徽省土壤黏粒含量大致分布為由北向南逐漸降低。利用RF模型進行預測制圖基本上能夠反映大尺度區(qū)域上土壤屬性的空間分布。

表5 土壤屬性的RF建模精度評價

圖2 RF模型預測散點圖

圖3 安徽省土壤屬性實測值及預測值空間分布

3 結論

1)安徽省內,對于土壤有機碳含量、土壤容重,土壤類型均是主要影響因子之一,可能是由于安徽省土地利用大部分為耕地,自然用地較少,導致人為因素對土壤屬性影響較大;對于土壤黏粒含量,高程和年均降水量為最主要的影響因素。

2)RF模型的建模結果表明,不同環(huán)境變量的組合分別解釋了研究區(qū)域內土壤有機碳含量、容重和黏粒含量的26%、23% 和22%,建模集和驗證集2相近,說明在大尺度區(qū)域內,RF模型能夠有效地減少過擬合問題且對于土壤屬性空間分布的預測具有較高的穩(wěn)定性。

3)土壤容重和黏粒含量的預測精度不是很高,原因可能是由于研究區(qū)域面積過大,不同區(qū)域地貌和氣候差異較大,以及一些可能影響土壤屬性的環(huán)境變量并沒有考慮到模型中。在以后的研究中可以增加采集樣本數(shù)量并加入更多的環(huán)境因子作為預測變量以提高預測精度。

[1] Lal R, Kimble J M, Stewart B A, et al. Global climate change and pedogenic carbonate[J]. Geoderma, 1999, 104(1): 135–141

[2] Dixon J B. Roles of clays in soils[J]. Applied Clay Science, 1991, 5(5/6): 489–503

[3] Kuhn M. Building predictive models in R using the caret Package[J]. Journal of Statistical Software, 2008, 28(5): 1–26

[4] Mansuy N, Thiffault E, Paré D, et al. Digital mapping of soil properties in Canadian managed forests at 250 m of resolution using the-nearest neighbor method[J]. Geoderma, 2014, s 235/236(4): 59–73

[5] Henderson B L, Bui E N, Moran C J, et al. Australia-wide predictions of soil properties using decision trees[J]. Geoderma, 2005,124(3): 383–398

[6] Breiman L. Random forests[J]. Machine learning, 2001, 45(1): 5–32

[7] Rodriguez-Galiano V F, Chica-Rivas M. Evaluation of different machine learning methods for land cover mapping of a Mediterranean area using multi-seasonal Landsat images and Digital Terrain Models[J]. International Journal of Digital Earth, 2014, 7(6): 492–509

[8] Dharumarajan S, Hegde R, Singh S K, et al. Spatial prediction of major soil properties using Random Forest techniques—A case study in semi-arid tropics of South India[J]. Geoderma Regional, 2017, 10: 154–162

[9] Chagas C D S, Junior W D C, Bhering S B, et al. Spatial prediction of soil surface texture in a semiarid region using random forest and multiple linear regressions[J]. Catena, 2016, 139: 232–240

[10] 郭澎濤, 李茂芬, 羅微, 等. 基于多源環(huán)境變量和隨機森林的橡膠園土壤全氮含量預測[J]. 農業(yè)工程學報, 2015, 31(5): 194–202

[11] 姜賽平, 張懷志, 張認連, 等. 基于三種空間預測模型的海南島土壤有機質空間分布研究[J]. 土壤學報, 2018, 55(4): 1007–1017

[12] 李德成, 張甘霖, 王華, 等. 中國土系志·安徽卷[M]. 北京: 科學出版社, 2017: 3–24

[13] 李欣海. 隨機森林模型在分類與回歸分析中的應用[J]. 應用昆蟲學報, 2013, 50(4): 1190–1197

[14] 李龍, 姚云峰, 秦富倉, 等. 基于地理加權回歸模型的土壤有機碳密度影響因子分析[J]. 科技導報, 2016, 34(2): 247–254

[15] 趙明松, 張甘霖, 李德成, 等. 江蘇省土壤有機質變異及其主要影響因素[J]. 生態(tài)學報, 2013, 33(16): 5058– 5066

[16] Chai T, Draxler R R. Root mean square error (RMSE) or mean absolute error (MAE)? - Arguments against avoiding RMSE in the literature[J]. Geoscientific Model Develop-ment Discussions, 2014, 7(3): 1247–1250

[17] Miller F P, Vandome A F, Mcbrewster J. Coefficient of determination[J]. Alphascript Publishing, 2006, 31(1): 63–64

[18] Gelman A. Analysis of variance[J]. Quality control & applied statistics, 2006, 20(1): 295–300

[19] Sonobe R, Tani H, Shimamura H, et al. Parameter tuning in the support vector machine and random forest and their performances in cross-and same-year crop classification using TerraSAR-X[J]. International Journal of Remote Sensing, 2014, 35(23): 7898–7909

[20] 許信旺, 潘根興, 曹志紅, 等. 安徽省土壤有機碳空間差異及影響因素[J]. 地理研究, 2007, 26(6): 1077–1086

[21] 趙明松, 李德成, 王世航. 近30年安徽省耕地土壤有機碳變化及影響因素[J]. 土壤學報, 2018, 55(3): 595–605

Spatial Prediction of Soil Properties Based on Random Forest Model in Anhui Province

LU Hongliang1, ZHAO Mingsong1,2*, LIU Binyin1, ZHANG Ping1, LU Longmei1

(1 School of Geodesy and Geomatics, Anhui University of Science and Technology, Huainan, Anhui 232001, China; 2 State Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008, China)

It is important to study the spatial variability and distribution of soil properties for understanding ecosystems, formulating agricultural policies, conducting soil management and monitoring environmental changes caused by land use. The purpose of this paper is to explore the accuracy of the spatial prediction of soil properties at the provincial scale by the Random Forest (RF) model. Anhui Province in East China was selected as the study area, soil data obtained during the 2ndNational Soil Survey and during 2010—2011 were used, the environmental variables were collected with GIS spatial analysis technique, and the correlation between environmental factors and soil properties was analyzed by RF model. The results showed thatin the RF modeling process, SOC prediction model was the most robust and the prediction accuracy was the highest when the mtry value was 1 and the ntree value was 1 000; when the mtry value was 1 and the ntree value was 1 000 and 100 respectively, soil bulk density (BD) and clay content prediction models were the best. The elevation, NDVI, landform, muti-resolution index of valley bottom flatness (MrVBF) and soil type were the most important predictors of SOC content; Landform, mean annual precipitation (MAP), MrVBF, elevation and soil type were the most important prediction factors of soil BD; Elevation, MAP, MrVBF and plan curvature were the most important predictors of soil clay content; RF model can be used for spatial prediction of soil properties and has certain advantages in treating the qualitative variables such as soil type and landform; Multi-source environmental variable combinations explained 26% of SOC content, 23% of soil Bd and 22% of clay content, respectively.The use of machine learning for predicting soil properties and digital soil mapping is more efficient than traditional methods, it is of significance to use RF model in spatially predicting soil properties in the large-scale area.

Soil properties prediction; Random Forest model; Environmental variables; Anhui Province

國家自然科學基金項目(41501226)、安徽省高校自然科學研究項目(KJ2015A034)和土壤與農業(yè)可持續(xù)發(fā)展國家重點實驗室開發(fā)基金項目(Y412201431)資助。

(zhaomingsonggis@163.com)

盧宏亮(1993—),男,安徽銅陵人,碩士研究生,主要從事數(shù)字土壤制圖研究。E-mail: 15656232332@163.com

S159

A

10.13758/j.cnki.tr.2019.03.025

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 色香蕉影院| 伊人久久影视| 国产色爱av资源综合区| 亚洲视频在线观看免费视频| 国产亚洲精品无码专| 日本一区二区三区精品国产| 久久天天躁夜夜躁狠狠| 国产91无毒不卡在线观看| 亚洲香蕉久久| 久久精品一品道久久精品| 制服丝袜亚洲| 99这里只有精品6| 久久这里只精品国产99热8| 亚洲一级无毛片无码在线免费视频| 国产乱人伦AV在线A| 国产尹人香蕉综合在线电影| 3344在线观看无码| 在线播放国产99re| 午夜欧美理论2019理论| 亚洲精品麻豆| 国产99欧美精品久久精品久久| 天天色综网| 国产成人无码综合亚洲日韩不卡| vvvv98国产成人综合青青| 91精品专区| 免费无码一区二区| 久久香蕉国产线看观看亚洲片| 欧美中出一区二区| 久久青草热| 精品三级在线| 国产视频你懂得| 蜜臀AV在线播放| 亚洲欧洲美色一区二区三区| 国产精品私拍在线爆乳| 精品日韩亚洲欧美高清a| 经典三级久久| 亚洲免费三区| 亚洲午夜福利在线| 色噜噜综合网| 一本大道香蕉高清久久| 啪啪永久免费av| 久久精品国产免费观看频道| 综合色区亚洲熟妇在线| 97在线碰| 午夜在线不卡| 69av在线| 中文字幕佐山爱一区二区免费| 欧美人人干| 亚洲人成成无码网WWW| 91亚洲精品第一| 国产精品美女自慰喷水| 日日噜噜夜夜狠狠视频| 国产玖玖视频| 国产精品尤物在线| 91视频区| 亚洲第一中文字幕| 国产在线拍偷自揄观看视频网站| 在线精品视频成人网| 国产在线拍偷自揄观看视频网站| 国产人碰人摸人爱免费视频| 成人毛片免费在线观看| 精品视频在线一区| 蝌蚪国产精品视频第一页| 欧美人与动牲交a欧美精品| aa级毛片毛片免费观看久| 亚洲精品波多野结衣| 这里只有精品在线播放| 91香蕉视频下载网站| 亚洲热线99精品视频| 中文字幕 日韩 欧美| 91精品啪在线观看国产91| 亚洲男人的天堂久久香蕉| 无码专区国产精品第一页| 91网在线| 依依成人精品无v国产| 国产精品一区二区在线播放| 亚洲高清无码精品| 日韩高清一区 | 国产成人综合久久精品尤物| 国产精彩视频在线观看| 天堂岛国av无码免费无禁网站| 亚洲色图欧美激情|