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

基于RF和SGT算法的子區(qū)優(yōu)先建模對綠洲尺度土壤鹽度預(yù)測精度的影響

2018-12-27 01:59:38王飛楊勝天魏陽楊曉東丁建麗1
中國農(nóng)業(yè)科學(xué) 2018年24期
關(guān)鍵詞:模型

王飛,楊勝天,魏陽,楊曉東,3,丁建麗1,,3

?

基于RF和SGT算法的子區(qū)優(yōu)先建模對綠洲尺度土壤鹽度預(yù)測精度的影響

王飛1,2,3,楊勝天2,魏陽2,楊曉東2,3,丁建麗1,2,3

(1新疆大學(xué)智慧城市與環(huán)境建模新疆普通高校重點實驗室,烏魯木齊 830046;2新疆大學(xué)資源與環(huán)境科學(xué)學(xué)院,烏魯木齊 830046;3綠洲生態(tài)教育部重點實驗室,烏魯木齊 830046)

【目的】試圖通過優(yōu)先在干旱區(qū)綠洲的子區(qū)構(gòu)建模型以提高綠洲全局土壤鹽度的預(yù)測精度。同時量化全局模型和子區(qū)模型之間精度的差異性和不確定性。【方法】利用隨機森林(Random Forest,RF)和隨機梯度增進算法(Stochastic Gradient Treeboost,SGT)定量化上述不確定性,同時,對比本地尺度多個情景(景觀)優(yōu)先建立模型再合并預(yù)測值對于模擬全局土壤鹽度的精度影響。基于驅(qū)動因子(土地利用和地貌),響應(yīng)因子(Normalized Difference Vegetation Index, NDVI和土壤電導(dǎo)率,EC),研究設(shè)計了27個能夠相對覆蓋典型綠洲不同土壤鹽度變異性的環(huán)境情景。【結(jié)果】70.37%(19/27)的情景證明SGT的預(yù)測精度高于RF。單獨建模的10個情景的預(yù)測精度高于全局模型下10個再分類情景(根據(jù)情景設(shè)定規(guī)則將全局模型預(yù)測值再分類)的精度。特別是,EC≤4 ds·m-1和 2 dS·m-1< EC<16 dS·m-1兩個情景應(yīng)該單獨進行建模預(yù)測。4個情景(兩兩合并)預(yù)測值合并后的精度高于全局模型再分類后的精度。需要指出的是,用于綠洲尺度子區(qū)情景構(gòu)建的首選分割變量是EC,其次是地貌和土地利用。【結(jié)論】研究推薦基于SGT在綠洲內(nèi)部不同景觀尺度上優(yōu)先建模,再將各景觀尺度的預(yù)測值進行合并,以提高綠洲土壤鹽度的推理精度。

土壤鹽分;機器學(xué)習(xí);干旱區(qū);Landsat OLI;空間異質(zhì)性;隨機森林算法;隨機梯度增進算法

0 引言

【研究意義】土壤鹽漬化的精準(zhǔn)預(yù)測是評估綠洲生態(tài)環(huán)境的適應(yīng)性和提高資源利用效率的重要前提。【前人研究進展】已有諸多基于遙感信息識別不同區(qū)域土壤鹽漬化影響程度和分布范圍的研究[1-8]。以上研究基于不同傳感器(IKONOS,Landsat, MODIS)獲取的數(shù)據(jù),計算與土壤鹽度相關(guān)的指數(shù),并建立全局土壤鹽度模型,成功預(yù)測了各自研究區(qū)內(nèi)鹽漬化土壤的分布范圍。【本研究切入點】然而,綠洲尺度下與土壤鹽分呈現(xiàn)響應(yīng)關(guān)系的環(huán)境變量是否可以在綠洲景觀組分內(nèi)(如土利用類型、地貌類型、不同植被覆蓋度)持有同等顯著性的空間變異解釋力尚不明確。綠洲內(nèi)影響陸氣水熱交換的因素(如地表粗糙度、反照率)隨著自然生境向人工灌溉農(nóng)業(yè)的轉(zhuǎn)變發(fā)生改變[9],致使不同土地利用類型之間表面能量和凈輻射量發(fā)生變化[10]。綠洲內(nèi)環(huán)境建構(gòu)組分如不同土地利用類型(農(nóng)田、灌木、草地、鹽堿地等)對水鹽運移的驅(qū)動機制也各有不同[11-13]。若僅依靠全局尺度建構(gòu)的土壤鹽漬化-環(huán)境變量模型模擬組分內(nèi)的土壤鹽度空間變異性,可能存在值域被夸大或被低估,甚至被誤判的現(xiàn)象(基于流域/綠洲尺度選擇的環(huán)境變量可能無法充分解釋組分內(nèi)部相對均勻環(huán)境中的土壤鹽度變化)[14-15]。組分尺度(農(nóng)田、灌木、草地、鹽堿地等)上對于土壤鹽漬化-環(huán)境響應(yīng)關(guān)系精確定量描述可能有助于提高綠洲尺度的土壤鹽漬化制圖精度。統(tǒng)計學(xué)方法是土壤鹽漬化定量研究中的重要手段。機器學(xué)習(xí)被定義為基于現(xiàn)代統(tǒng)計方法挖掘目標(biāo)和響應(yīng)變量之間關(guān)系模式的過程[16]。通過查找文獻,隨機森林(Random Forest, RF)和隨機梯度增進算法(Stochastic Gradient Treeboost, SGT)是目前土壤制圖領(lǐng)域應(yīng)用較為廣泛的機器學(xué)習(xí)方法[17-20]。RF的優(yōu)勢在于具備非線性挖掘能力;數(shù)據(jù)分布不需要符合任何假設(shè);同時處理類型和連續(xù)變量;防止過度擬合;有效抑制數(shù)據(jù)中存在的噪聲;定量描述變量的貢獻度;只需要率定少量參數(shù)[21]。SGT由CART改進而來,通常可以實現(xiàn)準(zhǔn)確和穩(wěn)健的預(yù)測,有效抑制過度擬合效應(yīng)[22]。不需要變量的先驗假設(shè),比傳統(tǒng)的廣義線性或加權(quán)模型提供更大的靈活性。存在空間異質(zhì)性和異常值時,SGT依然能獲得較高的預(yù)測精度。在復(fù)雜關(guān)系定義的干旱區(qū)土壤生態(tài)系統(tǒng),上述算法具備的優(yōu)勢顯得尤為重要。然而,這兩種優(yōu)秀的現(xiàn)代統(tǒng)計學(xué)算法鮮有用于預(yù)測干旱區(qū)土壤鹽度的先例。【擬解決的關(guān)鍵問題】本研究一是比較不同環(huán)境建構(gòu)組分下RF和SGT的預(yù)測能力;二是定量化描述綠洲分區(qū)建模對綠洲尺度鹽度預(yù)測精度的影響;三是獲取不同環(huán)境建構(gòu)組分條件下的敏感因子及其貢獻度。

1 材料與方法

1.1 研究區(qū)

渭干河-庫車河流域(渭庫)是位于塔里木盆地北部的一個典型綠洲。該研究區(qū)的中心經(jīng)緯度坐標(biāo)為82.50°N,41.38°E,主要地貌類型有:低海拔沖積扇平原、低海拔沖積平原、中等高度沖積扇平原、低海拔固定草灌木區(qū),以及低海拔半固定草灌木區(qū)(圖1)。研究區(qū)海拔范圍為892—1 100 m,從西北到東南逐漸降低。該區(qū)主要的土壤類型為:積土成土、鹽堿土、石膏層/鈣鹽土、鈣積變性土、鈣成土和石灰性黑土。灌區(qū)外圍受土壤鹽漬化和土壤水分的影響,整體的植被覆蓋度較低,自然植被主要包括:and[23]。該區(qū)為極端干旱荒漠氣候,年平均降雨量為51.6 mm,年平均蒸散量為2 356 mm,年積溫為(>10℃)為4 500℃。由于水庫引水灌溉急劇增加(由于灌溉土地的擴張),該區(qū)地下水位和地下水礦化度顯著增加。地下水位上升和地下水鹽化是造成土壤鹽漬化的主要因素[24]。耕地鹽漬化面積超過50%,其中30%的耕地受到嚴(yán)重影響[25]。

圖1 渭干河-庫車河流域景觀地貌特征及綠洲采樣點分布圖

1.2 數(shù)據(jù)與方法

1.2.1 土壤樣品采集與土壤鹽分分析處理 為了進行土壤鹽度分析和建模預(yù)測,本研究共計獲得了371個樣本。根據(jù)實際情況最終樣點的布局采取兩種不同的設(shè)計方式,采樣點分布情況見圖1。研究初步的樣點布局是基于限制拉丁超立方體(Conditioned Latin Hypercube Method,cLHS)而設(shè)計[26]。該方法是在已確定采樣數(shù)目的前提下抽取盡可能全面地覆蓋由定性或定量變量構(gòu)成的屬性空間的樣本。cLHS的程序和算法可在Minasny和McBratney(2006)文獻中可以找到。然而,在實際采樣中,由于各種原因(道路不通、植被過密或地形崎嶇),3/5的位置未能被抵達。因此,我們隨后采用了分層隨機抽樣的策略,樣點的布局綜合考慮了綠洲景觀類型、水文特征、道路通達性和土壤鹽漬化空間異質(zhì)性、土地利用和地形等因素。數(shù)據(jù)收集時間為2016年9月9—30日。每個樣點用五點梅花的方式取土,采樣深度為0—10 cm,隨后將測試的數(shù)據(jù)進行平均,作為本樣點的實際觀測值。

本研究中,土壤鹽度(測定土水比1﹕5溶液電導(dǎo)率,EC1:5)在實驗室內(nèi)進行分析,參照《土壤和農(nóng)業(yè)化學(xué)分析》,將土樣風(fēng)干后磨碎,過2 mm篩子,制備土壤溶液,比例為土水比1﹕5,溫度設(shè)定為25℃,用以測定土壤樣本的電導(dǎo)率(EC1:5)[27]。所有樣點的電導(dǎo)率平均值作為采樣點的基本值。此外,本研究區(qū)EC1:5和土壤鹽分(g·kg-1)具有高度相關(guān)性(2= 0.95),所以EC1:5可用作評估研究區(qū)土壤鹽分的替代參數(shù)。

1.2.2 Landsat OLI圖像預(yù)處理 本研究所使用的衛(wèi)星影像是Landsat OLI(145/31),獲取時間為2016年9月18日。基于ENVI 5.3中的FLAASH(Fast line- of-slight Atmospheric Analysis of Spectral Hypercubes)模型糾正大氣帶來的影響。糾正后的反射率數(shù)據(jù),其數(shù)值范圍為0—1。該數(shù)據(jù)用于計算后續(xù)研究中使用的各類土壤或植被相關(guān)指數(shù)。

1.2.3 基于驅(qū)動和響應(yīng)因子的場景設(shè)計 本研究基于驅(qū)動因子和環(huán)境響應(yīng)因子,設(shè)定以下12個場景用以量化分區(qū)建模對于綠洲尺度土壤電導(dǎo)率預(yù)測精度的影響(圖2)。將綠洲尺度的全部樣本構(gòu)建的模型設(shè)置為情景1。基于植被覆蓋情況,將植被覆蓋程度分為全覆蓋(NDVI>0.22,情景2)和無植被覆蓋(NDVI≤0.22,情景3)。NDVI閾值(0.22)的設(shè)定參考WU 等[5]和 SONG 等[28]的研究成果;基于土地利用:農(nóng)業(yè)用地(情景4)、草地(情景5)和未利用土地(情景6);基于地貌:低海拔的沖積扇平原(情景7,LAAFP)和低海拔固定灌木(情景8,LAFGS);基于土壤電導(dǎo)率分級(EC)[29]:EC ≤ 4 dS·m-1(情景9),2 dS·m-1<EC<16 dS·m-1(情景10),EC>4 dS·m-1(情景11),EC>16 dS·m-1(情景12)。場景的設(shè)置同時考慮了土壤發(fā)生學(xué)和人類活動影響。此外,研究還想說明分層預(yù)測非均質(zhì)地表土壤電導(dǎo)率的重要性。圖3顯示了每個場景的土地利用類型的組成。

圖2 綠洲-沙漠水熱和能量交換示意及情景覆蓋范圍(修改自LI等[10])

農(nóng)田:Farm land (FL);灌木 Shrub Lang(SL);高草地:High Grass Coverage (HGC);濕地:Wetland;低草地:Low Grass Coverage (LGC);沙地:Sandy Land (SYD);居民區(qū):Residential Area (RA);鹽堿地:Salt-affected Land (SAL);中草地:Medium Grass Coverage (MGC);果園:Orchard;戈壁:Gobi;裸巖地:Bare Rock Land (BRL)

1.2.4 環(huán)境變量 基于土壤方程選擇用于預(yù)測土壤電導(dǎo)率的環(huán)境變量:氣候因子(地表溫度)、生物因子(植被指數(shù))、DEM(數(shù)字高程模型)衍生變量集(30、60、90、120、150、180、210和240 m等8個分辨率)、母質(zhì)(地貌類型)、人類活動(土地利用)、土壤相關(guān)因子(波段、土壤相關(guān)指數(shù)、土壤濕度指數(shù))。具體變量詳見表1。

1.2.5 RF & SGT RF模型是基于決策樹發(fā)展而來的一種集成學(xué)習(xí)方法,通過多次 bootstrap抽樣獲得多個隨機樣本,并通過這些樣本分別建立相對應(yīng)的決策樹,從而構(gòu)成隨機森林。用于回歸的 RF,取所有決策樹預(yù)測結(jié)果的均值作為最終的預(yù)測結(jié)果[21]。RF算法主要涉及兩個關(guān)鍵參數(shù):一是對于每個參與模型的分支樹,即模型運算中的分裂次數(shù)mtry的設(shè)定。參與預(yù)測的變量個數(shù)不同,mtry的設(shè)定也不同,需要根據(jù)實際情況調(diào)整。本文針對設(shè)置的12個情景分別定義mtry。研究根據(jù)RF模型預(yù)測過程中(63%的樣本用于建模)產(chǎn)生的袋外誤差(37% 的樣本用于驗證)對mtry的值再做具體的調(diào)整。每個情景的mtry值由1到30分別進行測試,取袋外誤差最小的mtry值作為本情景最適宜的分裂次數(shù)。另一個是樹模型運算中每次生成的樹的數(shù)量,ntree,模型的計算量與 ntree的值成正比,在 ntree增加并不能顯著提高模型預(yù)測能力的情況下,ntree的設(shè)定要盡可能小(一般>500),本文設(shè)定為1 000。RF模型的預(yù)測通過R語言的random forest工具包(randomForest)實現(xiàn)[35]。

表1 基于Landsat OLI(30 m)和DEM(8個空間分辨率)衍生的環(huán)境變量

SGT是回歸樹和Boosting的集成。Boosting的核心思想是:初始狀態(tài)下為每個訓(xùn)練樣本賦予一樣的權(quán)重值,每次迭代訓(xùn)練提高錯分樣本的權(quán)重,降低分對樣本的權(quán)重。迭代N次之后,得到N個弱分類器,最終通過權(quán)重加和的方式集合成強分類器[36]。SGT是Boosting算法框架的分支,該算法為了減少上次結(jié)果的殘差,在減少殘差梯度方向上不斷建立新的模型,直到誤差不在降低為止。該算法需要設(shè)置以下四個參數(shù):學(xué)習(xí)速率,抽樣比例,樹的最大化子節(jié)點個數(shù)和每次生成的樹的數(shù)量。學(xué)習(xí)速度越低,迭代次數(shù)越多,一般小于0.1,此研究中設(shè)置為0.01[36]。抽樣比例為1時,每次迭代的樣本集相同,小于1時,抽取的訓(xùn)練樣本集都不同,有助于抑制過擬合,取值范圍0.5—0.8[17],為了和RF進行對比,這里設(shè)置為0.63。樹的最大化子節(jié)點個數(shù)需要根據(jù)實際情況設(shè)定,測試值的范圍1—15,取RMSE最小的值為最佳最大化子節(jié)點個數(shù)(圖6)。每次生成的樹的數(shù)量設(shè)置為1000。

1.2.6 變量篩選與最佳組合的選擇 研究表明通過去除潛在不相關(guān)的環(huán)境變量會提升模型的魯棒性。原因在于篩選過程可不斷降低不相關(guān)性的變量對模型精度的負(fù)向擾動,以提高預(yù)測的準(zhǔn)確性。本研究會在以上12種情景下分別進行最佳變量組合的測試。篩選步驟主要參考SVETNIK等[37]和HEUNG等[18]的研究成果。(1)基于RF和SGT,對變量的重要性進行排序。(2)依次刪除排在末尾最不重要的變量,并再進行迭代運算,驗證采用5-fold交叉驗證獲取RMSE。根據(jù)RMSE最低值確定最佳變量組合。(3)在各個設(shè)定的情景中重復(fù)步驟1和步驟2,確定每個場景的最優(yōu)變量組合。

1.2.7 情景對比方案 研究采用以下兩種方法比較情景(子區(qū))重組對綠洲尺度鹽度預(yù)測的影響。方案一基于全樣本構(gòu)建綠洲尺度預(yù)測模型,之后將預(yù)測值依據(jù)12個情景的設(shè)定規(guī)則重新組合(Regroup Scenarios,RS)(圖4-b)。另外,為了對比,如圖4-a所示,在12個情景中單獨建立預(yù)測模型(Independent Scenarios,IS)。將上述兩種模式生成的情景驗證結(jié)果進行對比。方案二則組合不同場景的預(yù)測值,即NDVI>0.22 & NDVI≤0.22,EC≤4 dS·m-1& EC>4 dS·m-1,農(nóng)用地,草地&未使用的土地, LAAFP & LAFGS,然后比較上述組合方案(Combined Scenarios,CS)(圖5-a)與全樣本預(yù)測值根據(jù)上述組合方式形成對應(yīng)組合(Reclassified Scenarios from Oasis Scale Model,RSOSM )的精度差異性(圖5-b)。

圖4 獨立情景建模與全樣本建模下預(yù)測值根據(jù)情景重分類之間的精度對比方案

圖5 情景預(yù)測值的合并對于綠洲尺度土壤鹽度預(yù)測影響的對比方案

1.3 模型評估

研究基于5折交叉驗證(5-fold cross validation)獲取各情景下最優(yōu)模型并進行驗證。此驗證測試過程相對于簡單的訓(xùn)練-驗證比例而言需要花費更多的計算工作量,并適用于較小的數(shù)據(jù)集,結(jié)果可靠且具備無偏性[17]。訓(xùn)練數(shù)據(jù)集被隨機分成5個子集, 即4/5 的訓(xùn)練數(shù)據(jù)被用于模型訓(xùn)練,剩下的1/5用于模型驗證。每個抽樣過程重復(fù)20次,取其平均值。圖6和圖7顯示了每個情景下經(jīng)過4/5訓(xùn)練數(shù)據(jù)獲取的RF算法中mtry的最優(yōu)值(黑色放框)和 SGT算法中的最大化子節(jié)點個數(shù)的最優(yōu)值(黑色放框)。研究采用相關(guān)系數(shù)(2)、均方根誤差(Root Mean Square Error,RMSE)和預(yù)測偏差比(Residual Predictive Deviation,RPD)作為誤差評價指標(biāo)。RMSE 越接近零,RPD 值越大,預(yù)測精度越高。

2 結(jié)果

2.1 土壤鹽度變異性

表2顯示了不同情況下土壤鹽分的基本統(tǒng)計信息。情景1的統(tǒng)計分析表明, 該區(qū)土壤電導(dǎo)率最大值為 184.5 dS·m-1,最小為 0.14 dS·m-1,平均為 31.32 dS·m-1。根據(jù)美國鹽度實驗室的鹽度分類標(biāo)準(zhǔn)(0—2 dS·m-1為非鹽漬化;2—4 dS·m-1為輕度鹽漬化;4—8 dS·m-1為中度鹽漬化;8—16 dS·m-1為重度鹽漬化;>16 dS·m-1為極端鹽漬化),約67% 的樣品(EC>4 dS·m-1)受到土壤鹽漬化的影響。整體而言,不同情景的變異系數(shù)(coefficient of variation,CV)表明,該區(qū)0—10 cm土層的土壤鹽度呈現(xiàn)中度或高度空間變異性。(CV<0.1為低變異性,0.1<CV<1.0為中度變異性,CV>1.0為高度變異性)。

圖6 不同情景RF算法中mtry的最優(yōu)值

圖7 不同情景下SGT算法中最大化節(jié)點(node)的最優(yōu)值

表2 12個情景的土壤電導(dǎo)率(dS·m-1)統(tǒng)計特征

2.2 RF和SGT在土壤鹽度預(yù)測中的表現(xiàn)

表3和表4分別顯示了RS和IS模式下基于RF和SGT計算獲取的誤差水平。根據(jù)2,RMSE 和 RPD的統(tǒng)計結(jié)果:(1)RF和SGT在情景2、3和5中的預(yù)測水平相似。(2)70%的情景中(情景2、3、5、6、7、8、11和12),SGT的預(yù)測精度高于RF。2依次增加了6.06%、0、5.13%、23.40%、-4.00%、51.52%、18.75% 和 16.67%。RMSE依次降低了2.21%、0.005%、1.47%、11.73%、16.14%、12.49%、3.79%和3.20%。RPD值依次增加了0.44%、0、1.56%、13.67%、18.80%、13.82%、4.10%、3.33%。RF在情景9和10的預(yù)測精度高于SGT,2值提高了26.09%和34.21%。RMSE的值降低了4.82%和11.18%。RPD值增加了5.26%和12.50%。(3)所有情景的預(yù)測精度水平<0.01。表4的結(jié)果顯示:(1)除了情景2,其余情景的RF和SGT的預(yù)測精度水平相似。情景2中,相對于RF,SGT的R2值增加了46.67%,RMSE降低了11.61%,RPD值增加了13.56%。(2)50%的情景顯示SGT的預(yù)測精度高于RF。全樣本模式下,SGT的預(yù)測精度高于RF。另外,由于農(nóng)田樣本的高變異性,RF和SGT都無法完成預(yù)測。

表3 獨立情景下(IS模式)基于隨機森林和隨機增進算法的精度驗證

**:<0.01

表4 RS模式基于隨機森林和隨機增進算法的精度驗證(根據(jù)情景劃定規(guī)則重分類)

ns:沒有顯著性No significance

表5顯示了1.2.7節(jié)中方案二的對比結(jié)果。4個合并的情景顯示SGT的預(yù)測精度高于RF,情景2 和情景3的RF預(yù)測精度高于SGT。相比于全樣本(綠洲尺度)預(yù)測值的再分配情景(RSOSM)而言,上述情景預(yù)測值合并后,后者較前者而言,2值提升的最大值為0.08,最小提升了0.01。RMSE值降低的幅度最大值為0.08,最小值為0.01。RPD值提升的最大值為0.11,最小值為0.01。

與觀測值相比,RF和SGT的預(yù)測值都存在低值區(qū)高估和高值區(qū)低估的現(xiàn)象,導(dǎo)致殘差出現(xiàn)線性趨勢(圖8和圖9)。上述現(xiàn)象出現(xiàn)在各情景中。被高估的值出現(xiàn)于有人類管理活動的綠洲內(nèi)部地區(qū)(用低含鹽水灌溉),而被低估的值分布于綠洲的外圍(鹽分來自于上游灌區(qū)),即灌溉區(qū)鹽分在綠洲-荒漠交錯帶的聚集區(qū),以及渭庫綠洲的東北角。

表5 情景合并模式(CS模式與RSOSM模式)對土壤鹽度預(yù)測的精度影響

**:<0.01

圖8 不同情景下基于RF預(yù)測的殘差與觀測值的對比

圖9 不同情景下基于SGT預(yù)測的殘差與觀測值的對比

2.3 情景分區(qū)建模的預(yù)測精度

表3和表4顯示了1.2.7節(jié)中方案一獲得的統(tǒng)計結(jié)果:(1)對于RF來說, 以下6個情景IS模式(表 3)的預(yù)測精確高于RS模式(表 4):情景2,5,6,9,10,12。以 RPD 值為例,上述方案的精度分別提升了4.23%、3.23%、2.21%、1900%、620% 和11.11%。其余4個情景的RS模式(表3)的預(yù)測精確高于IS模式:情景3、11、7、8。RPD 值增加了3.14%、2.4%、13.97% 和2.38%。(2)對于SGT來說,以下6個情景IS模式(表3)的預(yù)測精確高于RS模式(表4):情景5、6、8、9、10、12。RPD 值分別增加了4.84%、8.22%、2180%、540%、12.73% 和12.00%。反之,RS模式(表3)的預(yù)測精確高于IS模式的情景為:情景2、3、7。RPD 值增加5.97%、2.38% 和2.11%。

表5顯示了2.2.7節(jié)中方案二中的統(tǒng)計結(jié)果(圖5)。首先,將情景9和11的預(yù)測值合并,CS與RSOSM模式結(jié)果對比:基于SGT,2和 RPD 分別提升了14.59% 和6.47%,基于RF,分別提升了10.87% 和4.38%。其次,當(dāng)情景2和3的預(yù)測值合并后,CS與RSOSM模式結(jié)果相比:基于SGT,2和 RPD 分別降低了2.17%和2.90%,基于RF,二者分別降低了2.08%和1.46%。地貌分組當(dāng)中的情景7和8合并后,CS與RSOSM模式結(jié)果相比:SGT的預(yù)測結(jié)果顯示,前者模式中2值上升了8.51%,RPD下降了3.62%,RF預(yù)測結(jié)果顯示,2下降了 7.00%,RPD 值下降了3.03%。當(dāng)草地和未用土地的預(yù)測值組合時,CS與RSOSM模式結(jié)果對比,對于SGT,前者相比后者,預(yù)測精度增加17.50%(2)和5.38%(RPD),基于RF,預(yù)測精度提高了7.50%(2)和2.32%(RPD)。

2.4 不同情景分區(qū)指示土壤鹽度變化的重要變量

圖10顯示了10個情景模式下各自最優(yōu)變量集。RF和SGT都可產(chǎn)生變量的貢獻度,這里只顯示精度相對較高的算法獲取的變量及其貢獻度。將變量的貢獻度進行歸一化處理后經(jīng)對比發(fā)現(xiàn),由于外部環(huán)境或者所處地理位置的影響,各情景的關(guān)鍵變量組有一定的差異性。同時,研究根據(jù)變量計算的數(shù)據(jù)來源,重新核算了影響鹽漬化空間差異性的貢獻度。結(jié)果顯示,基于地形衍生的變量組(多空間分辨率)的貢獻度最高,平均值為52.98%,其次為Landsat衍生變量組,平均值為36.22%,地貌類型的平均貢獻度為11.29%,土地利用的平均貢獻度為10.51%。對比所有情景發(fā)現(xiàn),僅從變量出現(xiàn)的頻次考慮,地表溫度,EVI2,EVI,ENDVI,F(xiàn)SEN依次出現(xiàn)了7次、6次、6次、4次、4次。另外,地貌和土地利用在下列情景中的貢獻度位居前列:情景1、2、3、11、12。

圖10 不同情景下基于1.2.6節(jié)中介紹的方法迭代獲取的重要變量

2.5 土壤鹽度的空間分布特征

根據(jù)第1.2.6節(jié)的方法,獲取RF和SGT兩種算法迭代產(chǎn)生的模型,并繪制了研究區(qū)土壤鹽漬化的空間分布圖(圖11)。圖11中的a和b圖顯示了全樣本模型的預(yù)測結(jié)果,c和d,e和f,g和h分別對應(yīng)情景7和8,情景3和2,情景5和6。地貌的2個類型與土地利用中2個類型各自合并后,并不能覆蓋全綠洲,合并后外圍的數(shù)據(jù)則由全樣本模型的預(yù)測值進行填充。整體而言,從定性層面來看,上述結(jié)果與DING等[23]與WANG等[25]結(jié)果相近。綠洲內(nèi)部,灌區(qū)農(nóng)業(yè)種植區(qū)的土壤鹽度相對較低,由于地形地勢的引導(dǎo),鹽分都積聚在綠洲-荒漠過渡帶地區(qū)。

3 討論

3.1 RF和SGT的精度評價

觀察表3、4和5后發(fā)現(xiàn),27個情景中的18個(接近70%)都顯示 SGT的預(yù)測精度高于RF。此結(jié)果暗示,SGT相對更適合預(yù)測干旱區(qū)土壤鹽度空間變異性。然而,上述兩種算法在鹽漬化研究中并無比較的先例,在此僅引用其他研究加以說明。NAGHIBI等基于SGT,CART和RF研究阿富汗地區(qū)地下水噴泉潛在分布區(qū),結(jié)果顯示,SGT的預(yù)測效果最佳,其次為CART和RF[38]。YOUSSEF等[39]基于廣義線性模型(Generalized Linear Models,GLM),CART,SGT和RF評價沙特阿拉伯地區(qū)滑坡災(zāi)害危險性,AUC(Area Under the Curve ) 曲線結(jié)果顯示,SGT的值最高,為0.958,其次為GLM,為0.821,CART為0.816,最后為RF,0.783。值越大代表精度越高。YANG等[40]比較SGT與RF預(yù)測青藏高原東北部高植被覆蓋區(qū)土壤有機質(zhì)的空間分布,結(jié)果顯示,SGT的預(yù)測精度稍高于RF,二者預(yù)測的結(jié)果在空間分布上較為相似。

a和b:全樣本模型;c和d是情景7和情景8合并的結(jié)果;e和f是情景2和情景3合并的結(jié)果;g和h是情景5和情景6合并的結(jié)果

基于上述兩種算法的模擬值出現(xiàn)低值高估和高值低估的現(xiàn)象已出現(xiàn)在其他研究中[41-42]。主要的原因在于小數(shù)據(jù)集樣本的代表性不足。本研究的樣本設(shè)計起初是以綠洲尺度為前提而定,并不是完全基于小尺度而定。聚焦至不同情景下,樣點的代表性可能有所下降,致使該模式(全樣本建模)下并不是每個情景的預(yù)測精度都可以達到綠洲尺度的水平。上述現(xiàn)象可能無法完全避免,在于研究區(qū)土壤鹽度空間變異性較強。樣本的設(shè)計依據(jù)來自于環(huán)境變量,這些變量在空間分布上并不能完全與土壤屬性呈高度一致性變化。對此,面對復(fù)雜的變異環(huán)境需要加入更多不同類型且更為普及的高空間分辨率傳感器進行相互校正或能提高樣本的代表性,以及滿足大尺度制圖的需求。

基于Landsat 數(shù)據(jù)建立土壤鹽度模型的案例已有報道(表1)。涉及的深度包括0—10、0—20和0—30 cm,部分研究未列出深度。參與建模的樣本量殘次不齊。選擇的變量或為單一變量或為變量組。相關(guān)性有高有低(2= 0.874,2= 0.564,2= 0.483,2= 0.78,2= 0.45,2= 0.71,2= 0.93,等)。選擇的方法各有優(yōu)缺點(線性、多元線性方程、神經(jīng)網(wǎng)絡(luò)、支持向量機、指數(shù)模型等)。驗證模型包括設(shè)定訓(xùn)練樣本/驗證樣本的比例和十字交叉驗證。同時,各自的地理環(huán)境也有明顯的差異性。由于上述多維信息的差異性,不能直接對比。但本研究認(rèn)為所得結(jié)果進一步豐富了該領(lǐng)域關(guān)于環(huán)境-土壤鹽漬化關(guān)系知識庫。

3.2 情景分區(qū)建模對綠洲土壤鹽度預(yù)測精度的影響

表3和表4中情景2、5、6、8、9、10、12,共計7個情景(占70%)的結(jié)果顯示了分區(qū)建模的有效性。其中EC≤4 dS·m-1和2 dS·m-1<EC<16 dS·m-1兩個情景的IS模式精度顯著高于RS模式。原因可能歸結(jié)于用于綠洲尺度建模的鹽度變異性與上述范圍的變異性相差太大,二者數(shù)據(jù)的值域范圍重疊度較低。此外,低鹽度區(qū)域如農(nóng)田,其土壤中鹽度的空間變異性與外界環(huán)境的變化同步性較其他地區(qū)并不明顯,控制該區(qū)植被生長的關(guān)鍵變量以土壤水分為主。同時,用于測試低鹽度的儀器在該區(qū)的測試精度并不穩(wěn)定。因此上述結(jié)果給我們的啟示在于綠洲尺度的模型并不適用于灌溉區(qū),需要在此情景下單獨建立相應(yīng)的預(yù)測模型。

表5中的組合場景再次驗證了樣本數(shù)據(jù)分割(根據(jù)地理特征)對提高土壤鹽度預(yù)測精度的重要性。場景9和場景11組合可以有效地提高預(yù)測精度,這是所有4種組合中精度提高最明顯的。對表5結(jié)果的綜合比較表明,當(dāng)有歷史土壤電導(dǎo)率數(shù)據(jù)時,它們是首選子區(qū)分割參考媒介。RF結(jié)果顯示了場景2和3被組合在一起時提高了預(yù)測精度。這表明植被覆蓋對土壤鹽度預(yù)測有一定的影響。然而,由于植被類型多樣,精度的提高范圍有限,因為一定程度的概率表明,同為植被覆蓋下有著相似NDVI值的單一或多種植被類型有著不同的土壤電導(dǎo)率值。利用地貌單元作為相對均勻區(qū)域的劃分基礎(chǔ)結(jié)果表明可以有效提高預(yù)測精度。這些結(jié)果的優(yōu)勢是分區(qū)范圍相對穩(wěn)定,極少隨時間變化。這與與母質(zhì)有關(guān),后者與土壤鹽漬化的發(fā)生和發(fā)展有間接關(guān)系[4,43]。土地利用數(shù)據(jù)的重要性(它整合了人類活動信息用以提高土壤鹽度預(yù)測的準(zhǔn)確性)排在情景9和11(參考2和RPD)的組合之后。這一發(fā)現(xiàn)也反映在各場景的重要變量排序中。值得注意的是,上述的基于地形和土地利用分區(qū)子區(qū)域,只覆蓋了整個綠洲的兩種類型,在下最終結(jié)論之前,我們還需要收集一定數(shù)量的樣本進一步驗證。

3.3 不確定性分析

全面分析上述結(jié)果后,未來的研究可能需要考慮以下幾個方面:(1)綠洲尺度下采樣點設(shè)計不僅要考慮整個研究區(qū),還要考慮子區(qū)域,并將道路可達性和采樣成本考慮在內(nèi)。雖然該研究在每個子區(qū)測試了模擬結(jié)果,預(yù)測誤差降低,如果在該區(qū)中收集更多代表性樣本,預(yù)測的準(zhǔn)確性將得到顯著改善。(2)由于混合像元問題,現(xiàn)存遙感或數(shù)字檔案數(shù)據(jù)與地理現(xiàn)象無法實現(xiàn)精確匹配。這個過程需要使用更高的空間分辨率數(shù)據(jù)。(3)在一定土壤鹽度值范圍內(nèi),植被物種多樣性和土壤的復(fù)雜性加劇了不確定性。(4)當(dāng)土壤中的鹽含量低于15%時,響應(yīng)變化不能從外界觀察得到,這就增加了兩種情況下的探查難度:一種是在農(nóng)業(yè)作物區(qū),另一種是深層地下水已經(jīng)停止積鹽。(5)遙感數(shù)據(jù),即使經(jīng)過了大氣校正和地形校正,在后續(xù)流程仍有一定的錯誤積累。

4 結(jié)論

在本研究中模擬的27個場景中(12個原始和15個衍生),70.37%的情景證明SGT比RF有更高的預(yù)測精度。因此,相對而言,在干旱區(qū)SGT比RF更適合于預(yù)測復(fù)雜環(huán)境下土壤鹽度的空間變異性。

將綠洲劃分為若干子區(qū),然后結(jié)合子區(qū)的預(yù)測值,能有效地提高綠洲尺度的預(yù)測精度。此結(jié)論特別適用于EC≤4 dS·m-1和2 dS·m-1<EC<16 dS·m-1兩個情景,需要單獨預(yù)測。首選分區(qū)媒介是土壤電導(dǎo)率(前提是數(shù)據(jù)能夠反映當(dāng)前土壤鹽度變異性),其次是地形和土地利用。

通過對多個情景優(yōu)化數(shù)據(jù)集的比較,以下變量對于指示子區(qū)或綠洲尺度的土壤鹽度預(yù)測具有重要的作用:土地利用、地貌、EVI、EVI2、FSEN、TEM、ENDVI以及多種空間分辨率的地形變量。

研究建議為了達到全組分景觀單元都保持較高的屬性預(yù)測精度,應(yīng)該在樣本設(shè)計上著手,并經(jīng)過反復(fù)測試找到最為合適的景觀分區(qū)變量,分區(qū)媒介的選擇可以從驅(qū)動或者響應(yīng)變量入手。

[1] ALLBED A, KUMAR L, ALDAKHEEL Y Y. Assessing soil salinity using soil salinity and vegetation indices derived from IKONOS high-spatial resolution imageries: Applications in a date palm dominated region.2014, 230:1-8.

[2] SCUDIERO E, SKAGGS T H, CORWIN D L. Regional scale soil salinity evaluation using Landsat 7, western San Joaquin Valley, California, USA.2014(2/3): 82-90.

[3] SCUDIERO E, SKAGGS T H, CORWIN D L. Regional-scale soil salinity assessment using Landsat ETM + canopy reflectance.2015, 169: 335-343.

[4] TAGHIZADEH-MEHRJARDI R, MINASNY B, SARMADIAN F, MALONE B P. Digital mapping of soil salinity in Ardakan region, central Iran.2014, 213: 15-28.

[5] WU W, MHAIMEED A S, AL-SHAFIE W M, ZIADAT F, DHEHIBI B, NANGIA V, PAUW E D. Mapping soil salinity changes using remote sensing in Central Iraq.2014(2/3): 21-31.

[6] YAHIAOUI I, DOUAOUI A, QIANG Z, ZIANE A. Soil salinity prediction in the Lower Cheliff plain(Algeria) based on remote sensing and topographic feature analysis.2015, 7(6): 794-805.

[7] ZHANG T T, QI J G, GAO Y, OUYANG Z T, ZENG S L, ZHAO B. Detecting soil salinity with MODIS time series VI data.2015, 52: 480-489.

[8] ZHANG T-T, ZENG S-L, GAO Y, OUYANG Z-T, LI B, FANG C-M, ZHAO B. Using hyperspectral vegetation indices as a proxy to monitor soil salinity.2011, 11(6): 1552-1562.

[9] GARCíA M, OYONARTE C, VILLAGARCíA L, CONTRERAS S, DOMINGO F, PUIGDEFáBREGAS J. Monitoring land degradation risk using ASTER data: The non-evaporative fraction as an indicator of ecosystem function.2008, 112(9): 3720-3736.

[10] LI X, YANG K, ZHOU Y. Progress in the study of oasis-desert interactions.2016, 230-231:1-7. DOI: 10.1016/j.agrformet.2016.08.022

[11] GONG L, RAN Q, HE G, TIYIP T. A soil quality assessment under different land use types in Keriya river basin, Southern Xinjiang, China.2015, 146: 223-229.

[12] TUTEJA N K, BEALE G, DAWES W, VAZE J. Predicting the effects of landuse change on water and salt balance-a case study of a catchment affected by dryland salinity in NSW, Australia.2003, 283(1): 67-90.

[13] WANG Y, LI Y. Land exploitation resulting in soil salinization in a desert–oasis ecotone.2013, 100: 50-56.

[14] HENGL T, MENDES D J J, HEUVELINK G B, RUIPEREZ G M, KILIBARDA M, BLAGOTI? A, SHANGGUAN W, WRIGHT M N, GENG X, BAUERMARSCHALLINGER B. SoilGrids250m: Global gridded soil information based on machine learning.2017, 12(2): e0169748.

[15] SCHILLACI C, LOMBARDO L, SAIA S, FANTAPPIè M, M?RKER M, ACUTIS M. Modelling the topsoil carbon stock of agricultural lands with the Stochastic Gradient Treeboost in a semi-arid Mediterranean region.2017, 286: 35-45.

[16] HASTIE T, TIBSHIRANI R, FRIEDMAN J H, FRANKLIN J. The elements of statistical learning, second edition: data mining, inference, and prediction.2009, 27(2): 83-85.

[17] ANGILERI S E, CONOSCENTI C, HOCHSCHILD V, M?RKER M, ROTIGLIANO E, AGNESI V. Water erosion susceptibility mapping by applying Stochastic Gradient Treeboost to the Imera Meridionale River Basin (Sicily, Italy).2016, 262: 61-76.

[18] HEUNG B, BULMER C E, SCHMIDT M G. Predictive soil parent material mapping at a regional-scale: A Random Forest approach.2014, 214: 141-154.

[19] HEUNG B, HO H C, ZHANG J, KNUDBY A, BULMER C E, SCHMIDT M G. An overview and comparison of machine-learning techniques for classification purposes in digital soil mapping.2016, 265: 62-77.

[20] LIE? M, GLASER B, HUWE B. Uncertainty in the spatial prediction of soil texture: Comparison of regression tree and Random Forest models.2012, 170: 70-79.

[21] BREIMAN L. Random Forests.2001, 45(1): 5-32.

[22] FRIEDMAN J H. Stochastic Gradient Boosting.2002, 38(4): 367-378.

[23] DING J, YU D. Monitoring and evaluating spatial variability of soil salinity in dry and wet seasons in the Werigan–Kuqa Oasis, China, using remote sensing and electromagnetic induction instruments.2014, 235-236: 316-322.

[24] BENNETT S J, BARRETTLENNARD E G, COLMER T D. Salinity and waterlogging as constraints to saltland pasture production: a review.2009, 129(4): 349-360.

[25] WANG F, CHEN X, LUO G, HAN Q. Mapping of regional soil salinities in xinjiang and strategies for amelioration and management.2015, 25(3): 321-336.

[26] MINASNY B, MCBRATNEY A B. A conditioned Latin hypercube method for sampling in the presence of ancillary information2006, 32(9): 1378-1388.

[27] 魯如坤. 土壤農(nóng)業(yè)化學(xué)分析方法. 北京: 中國農(nóng)業(yè)科技出版社, 1999.LU R K.Beijing: China Agricultural Science and Technology Press, 1999. (in Chinese).

[28] SONG W, MU X, RUAN G, GAO Z, LI L, YAN G. Estimating fractional vegetation cover and the vegetation index of bare soil and highly dense vegetation with a physically based method.2017, 58: 168-176.

[29] RICHARDS L A. Diagnosis and Improvement of Saline and Alkali Soils.1954, 60(3): 290.

[30] MONDAL P. Quantifying surface gradients with a 2-band Enhanced Vegetation Index (EVI2).2011, 11(3): 918-924.

[31] 陳紅艷, 趙庚星, 陳敬春, 王瑞燕, 高明秀. 基于改進植被指數(shù)的黃河口區(qū)鹽漬土鹽分遙感反演. 農(nóng)業(yè)工程學(xué)報, 2015, 31(5): 107-114.CHEN H Y, ZHAO G X, CHEN J C, WANG R Y, GAO M X. Remote sensing inversion of saline soil salinity based on modified vegetation index in estuary area of Yellow River.2015, 31(5): 107-114. (in Chinese)

[32] METTERNICHT G I, ZINCK J A. Remote sensing of soil salinity: potentials and constraints.2003, 85(1): 1-20.

[33] YU R, LIU T, XU Y, ZHU C, ZHANG Q, QU Z, LIU X, LI C. Analysis of salinization dynamics by remote sensing in Hetao Irrigation District of North China.2010, 97(12): 1952-1960.

[34] CECCATO P, GOBRON N, FLASSE S, PINTY B, TARANTOLA S. Designing a spectral index to estimate vegetation water content from remote sensing data: Part 1: Theoretical approach.2002, 82(2): 188-197.

[35] LIAW A, WIENER M. Classification and Regression by randomForest.2002(2/3): 18-22.

[36] ELITH J, LEATHWICK J R, HASTIE T. A working guide to boosted regression trees.2008, 77(4): 802-813.

[37] SVETNIK V, LIAW A, TONG C, CULBERSON J C, SHERIDAN R P, FEUSTON B P. Random Forest:? A Classification and Regression Tool for Compound Classification and QSAR Modeling.2003, 43(6): 1947-1958.

[38] NAGHIBI S A, POURGHASEMI H R, DIXON B. GIS-based groundwater potential mapping using boosted regression tree, classification and regression tree, and random forest machine learning models in Iran.2016, 188(1): 44.

[39] YOUSSEF A M, POURGHASEMI H R, POURTAGHI Z S, AL-KATHEERI M M. 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.2016, 13(5): 839-856.

[40] YANG R M, ZHANG G L, LIU F, LU Y Y, YANG F, YANG F, YANG M, ZHAO Y G, LI D C. Comparison of boosted regression tree and random forest models for mapping topsoil organic carbon concentration in an alpine ecosystem.2016, 60: 870-878.

[41] BLANCO C M G, GOMEZ V M B, CRESPO P, LIE? M. Spatial prediction of soil water retention in a Páramo landscape: Methodological insight into machine learning using random forest.2018, 316: 100-114.

[42] MULDER V L, LACOSTE M, RICHER-DE-FORGES A C, MARTIN M P, ARROUAYS D. National versus global modelling the 3D distribution of soil organic carbon in mainland France.2016, 263:16-34.

[43] 王玉剛, 李彥, 肖篤寧. 土地利用對天山北麓土壤鹽漬化的影響. 水土保持學(xué)報, 2009, 23(5): 179-183.WANG Y G, LI Y, XIAO D N. Effects of land use type on soil salinization at Northern Slope of Tianshan Mountain.2009, 23(5): 179-183. (in Chinese)

Influence of Sub-Region Priority Modeling Constructed by Random Forest and Stochastic Gradient Treeboost on the Accuracy of Soil Salinity Prediction in Oasis Scale

WANG Fei1,2,3, YANG ShengTian2, WEI Yang2, YANG XiaoDong2,3, DING JianLi1,2,3

(1Xinjiang Common University Key Laboratory of Smart City and Environmental Stimulation, Xinjiang University, Urumqi 830046;2College of Resource and Environmental Sciences, Xinjiang University, Urumqi 830046;3Laboratory for Oasis Ecosystem, Ministry of Education, Urumqi 830046)

【Objective】This study attempts to improve the prediction accuracy of soil salinity in arid oasis by building models preferentially in the sub-area of oasis. At the same time, the difference and uncertainty of accuracy between global model and sub-region model are quantified. 【Method】Therefore, to investigate the above differences, this study used two machine learning methods (Random Forest, RF and Stochastic Gradient Treeboost, SGT) to quantify the above effects and to prove the necessity of the building model in the sub-region compared with the full-sample model with respect to the simulation precision under the complex background of an arid region. Twenty-seven environmental scenarios (twelve original and fifteen derivatives) were designed based on the driving factors (land use and landform) and response factors (Normalized Difference Vegetation Index, NDVI and electrical conductivity, EC), which reflected variety of variabilities in soil salinity. After analyzing the results, the following preliminary conclusions were drawn. 【Result】The simulation results from 70.37% (19/27) of the scenarios showed that the predicted value of soil salinity from SGT was closer to the observed value from RF. Ten original sub-regions were modeled individually and compared with the full-sample model under the oasis scale (according to the 10 partition rules to reclassify the simulated values), and the result showed that the prediction accuracy of the former 70% scenario was higher than that of the latter. In particular, the regions of EC≤4 dS·m-1and 2 ddS·m-1<EC<16 dS·m-1should be modeled separately to predict the spatial variability of regional salinity. By combining the predictions of sub-regions and comparing them with the predicted values of the full-sample model, the former (all four different combination modes) showed a higher prediction accuracy than the latter. In addition, this result also indicated that the preferred medium for partitioning the sub-regions was soil electrical conductivity, followed by landform and land use. 【Conclusion】 The study proposes to establish a soil salinity model based on SGT preferentially on different landscape scales within the oasis, and then combine the predicted values of each landscape scale to improve the prediction accuracy of oasis soil salinity.

soil salinity; machine learning; arid regions; Landsat OLI; spatial heterogeneity; Random Forest; Stochastic Gradient Treeboost

2018-05-14;

2018-07-20

國家自然科學(xué)基金(U1603241、41661046、41771470、41261090、U1303381)、新疆大學(xué)博士研究基金(BS150246)

王飛,E-mail:volitation610@163.com。

丁建麗,E-mail:dingjianlixjdx@126.com

10.3864/j.issn.0578-1752.2018.24.007

(責(zé)任編輯 李云霞)

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产免费久久精品99re不卡| 日韩无码真实干出血视频| 国产精品黄色片| 亚洲欧美日韩精品专区| AV色爱天堂网| 波多野吉衣一区二区三区av| 久久免费精品琪琪| 91精品国产福利| 红杏AV在线无码| 性视频一区| 青青国产成人免费精品视频| 国产三级韩国三级理| 欧美日本视频在线观看| 欧美高清国产| 精品久久高清| 精品国产Av电影无码久久久| 日韩第九页| 亚洲综合一区国产精品| 老色鬼久久亚洲AV综合| 日韩av电影一区二区三区四区| Aⅴ无码专区在线观看| 欧美激情一区二区三区成人| 国产sm重味一区二区三区| 欧美色亚洲| 欧美在线精品怡红院| 女人18毛片水真多国产| 青草视频在线观看国产| 尤物在线观看乱码| 日本91在线| 欧美成人精品高清在线下载| A级毛片无码久久精品免费| 欧美日韩第三页| 欧美一级黄片一区2区| 91丝袜乱伦| 2020精品极品国产色在线观看| 国产尹人香蕉综合在线电影| 亚洲男人的天堂久久精品| 伊人91视频| 色悠久久久久久久综合网伊人| 久久人人爽人人爽人人片aV东京热 | 久久大香伊蕉在人线观看热2| 自拍亚洲欧美精品| 四虎永久免费地址在线网站| 国产精品第一区在线观看| 男女性色大片免费网站| 国产性猛交XXXX免费看| 国产99精品久久| 日本手机在线视频| a毛片免费观看| 乱色熟女综合一区二区| 综合久久五月天| 最新日韩AV网址在线观看| 国产自在线播放| 亚洲第一黄色网| 日韩在线播放中文字幕| 久久精品亚洲热综合一区二区| 久无码久无码av无码| 无码国产偷倩在线播放老年人| 亚洲成人免费在线| 午夜国产在线观看| 国产成人91精品| 婷婷99视频精品全部在线观看| 亚洲va视频| 97国产精品视频人人做人人爱| 国产00高中生在线播放| 午夜欧美在线| 日韩福利在线视频| 亚洲资源在线视频| 99热这里只有精品免费| 国产区精品高清在线观看| 在线欧美日韩| 国产SUV精品一区二区6| 婷婷久久综合九色综合88| 欧美精品xx| 国产SUV精品一区二区6| 欧美精品H在线播放| 97视频在线精品国自产拍| 青青草原国产一区二区| 欧美日韩一区二区在线免费观看| 亚洲无码高清免费视频亚洲| 国产黑丝一区| 欧美福利在线观看|