• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于地理加權(quán)回歸的地形平緩區(qū)土壤有機(jī)質(zhì)空間建模

    2019-12-19 01:13:18趙明松劉斌寅盧宏亮李德成張甘霖
    農(nóng)業(yè)工程學(xué)報 2019年20期
    關(guān)鍵詞:樣點(diǎn)回歸系數(shù)制圖

    趙明松,劉斌寅,盧宏亮,李德成,張甘霖

    基于地理加權(quán)回歸的地形平緩區(qū)土壤有機(jī)質(zhì)空間建模

    趙明松1,2,劉斌寅1,盧宏亮1,李德成2,張甘霖2※

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

    氣候變化效應(yīng)評估、土壤固碳潛力和肥力管理等,迫切需要詳盡的土壤有機(jī)質(zhì)(soil organic matter, SOM)空間分布信息。該文以江蘇省第二次土壤普查的1 519個典型土壤剖面的表層(0~20 cm)SOM含量為例,選擇1 217個樣本為建模集,302個為驗(yàn)證集,選取年均溫度、年均降雨、物理性黏粒和土壤pH值等因子進(jìn)行SOM的地理加權(quán)回歸(geographically weighted regression, GWR)建模。從建模集中分別隨機(jī)抽取100%(1 217個)、80%(973個)、60%(730個)、40%(486個),20%(243個)的樣點(diǎn),對比不同樣點(diǎn)數(shù)量下GWR和傳統(tǒng)全局回歸模型的精度差異,并選擇最優(yōu)模型進(jìn)行SOM空間預(yù)測制圖。結(jié)果表明:1)江蘇省SOM含量在不同空間尺度上存在極顯著的空間自相關(guān)性。不同樣點(diǎn)數(shù)量的建模集的全局自相關(guān)性和局部空間自相關(guān)聚類圖結(jié)果相似。全局Moran’s I值介于0.25~0.61(<0.001)。SOM含量空間分布以空間聚集特征為主,“高-高”聚集區(qū)主要分布在蘇中和蘇南地區(qū),“低-低”聚集區(qū)主要分布在蘇北地區(qū)。2)GWR建模結(jié)果均優(yōu)于傳統(tǒng)的全局回歸建模,其殘差在不同的空間尺度上均不存在空間自相關(guān)性。不同建模集的GWR的2adj較全局建模均提高0.15~0.20,其AIC和RSS均比全局模型有大幅降低,為56.08~360.19和17.40~76.67。不同建模樣本數(shù)量的GWR模型對SOM的解釋能力差異較小。3)建模樣點(diǎn)數(shù)量(除建模樣本=243)對GWR預(yù)測制圖結(jié)果的精度影響不大,RMSE介于5.56~5.75 g/kg之間,MAE介于3.87~4.05 g/kg之間,2介于0.52~0.48之間,均優(yōu)于全部建模樣點(diǎn)的普通克里格插值驗(yàn)證結(jié)果。該研究可為樣點(diǎn)數(shù)較少的省級尺度地區(qū)SOM空間建模與制圖提供借鑒。

    土壤; 有機(jī)質(zhì);模型;地理加權(quán)回歸;數(shù)字土壤制圖;地形平緩區(qū);江蘇省

    0 引 言

    土壤有機(jī)質(zhì)(soil organic matter,SOM)是土壤肥力和質(zhì)量的重要指標(biāo),是陸地生態(tài)系統(tǒng)碳儲量的重要組成部分,在農(nóng)業(yè)生產(chǎn)和環(huán)境保護(hù)等方面有著重要的意義[1-2]。研究SOM空間分布是評估區(qū)域土壤碳儲量、實(shí)現(xiàn)土壤可持續(xù)利用的前提,也是土壤質(zhì)量研究的重要內(nèi)容。20世紀(jì)90年代興起的計量土壤學(xué)和數(shù)字土壤制圖研究,主要基于土壤-景觀模型,利用土壤在地形、植被等景觀環(huán)境上的差異,采用線性回歸、廣義線性回歸、地統(tǒng)計學(xué)、模糊聚類、地理加權(quán)回歸(geographically weighted regression, GWR)、機(jī)器學(xué)習(xí)等數(shù)學(xué)方法和地理信息系統(tǒng)等技術(shù)進(jìn)行土壤屬性的空間建模,為獲取柵格化的土壤空間分布信息提供解決途徑[3-8]。建立較優(yōu)的土壤屬性空間預(yù)測模型是獲取較準(zhǔn)確的土壤空間分布信息的關(guān)鍵。但是在平原或地形平緩地區(qū),由于地形差異減小,植被等景觀環(huán)境因素差異也變小,地形和植被等條件不能反映土壤的空間變異,因此不能有效地用于土壤屬性的空間建模。

    土壤是具有高度變異的時空連續(xù)體,在形成發(fā)育過程中成土因子對土壤的作用是非線性,在較大區(qū)域中土壤屬性的主要影響因素較多,這種非線性作用更為明顯。周濤等[9]研究發(fā)現(xiàn)在<10 ℃的地區(qū),有機(jī)碳儲量與年均溫的負(fù)相關(guān)性最強(qiáng);在10~20 ℃的地區(qū),受與降水正相關(guān)的影響,有機(jī)碳儲量與年均溫呈正相關(guān);而在>20 ℃的地區(qū),有機(jī)碳儲量與溫度和降水的相關(guān)性都很差。趙明松等[10]研究表明江蘇省土壤質(zhì)地對SOM的作用在不同土壤類型和地理區(qū)域間存在著差異,總體上由北至南相關(guān)性減弱。因此在較大區(qū)域尺度上,利用傳統(tǒng)的全局模型(如線性回歸,普通回歸克里格等)并不能很好地或較為準(zhǔn)確地擬合環(huán)境因子與土壤屬性之間的非線性關(guān)系。

    GWR模型由英國地理學(xué)家Fotheringham等[11]提出,針對自變量和因變量間的關(guān)系在不同空間上的變化進(jìn)行非參數(shù)局部空間回歸建模,模型中自變量的回歸系數(shù)隨空間位置而變化,對于空間數(shù)據(jù)具有較強(qiáng)的局部分析能力[12-13]。GWR探測空間關(guān)系的非平穩(wěn)性,廣泛用于區(qū)域經(jīng)濟(jì)[14-16]、生態(tài)環(huán)境[17-20]、土壤建模制圖[21-25]等研究。劉瓊峰等[17]利用GWR模型分析了長沙城郊農(nóng)田土壤Pb、Cd含量的空間結(jié)構(gòu)及影響因素。郭龍等[26]等對比發(fā)現(xiàn)協(xié)同克里格和GWR模型在土壤預(yù)測制圖中均具有較高精度。Zeng等[25]研究表明在小流域尺度和地區(qū)尺度上GWR模型均能提高SOM建模的精度。Song等[21]、Yang等[23]研究表明GWR模型能夠提高較大復(fù)雜景觀區(qū)域的土壤有機(jī)碳和電導(dǎo)率的建模預(yù)測精度。Zhang等[24]基于GWR模型利用年均降雨、土壤類型和土地利用等因子預(yù)測了愛爾蘭全國的土壤有機(jī)碳的空間分布。上述研究表明了GWR模型在較大區(qū)域中土壤屬性空間建模的優(yōu)勢。

    針對較大區(qū)域尺度、地形平緩區(qū)土壤屬性空間建模中存在的問題,以及省級尺度上采集大量土壤樣品的困難。本研究選擇江蘇省為例,以第二次土壤普查數(shù)據(jù)為基礎(chǔ),選擇柵格氣象數(shù)據(jù)等,采用GWR模型開展地形平緩區(qū)SOM建模預(yù)測研究,并探討樣點(diǎn)數(shù)量對GWR建模精度的影響。研究結(jié)果可為省級尺度區(qū)域的SOM空間預(yù)測制圖提供借鑒,為江蘇省土壤肥力和質(zhì)量管理、提升區(qū)域土壤固碳能力等提供數(shù)據(jù)基礎(chǔ)。

    1 材料與方法

    1.1 研究區(qū)概況

    江蘇?。?16°18′~121°57′ E,30°45′~35°20′ N)地處中國大陸東部,長江、淮河下游,面積10.26萬km2。全省處于亞熱帶向暖溫帶過渡地帶,年均氣溫13~16 ℃,年均降水量800~1 200 mm,年均日照2 100~2 600 h。全省平均海拔13 m,海拔低于10 m或坡度小于1%的區(qū)域占全省陸地面積的90%。按地質(zhì)地貌差異,分為低山丘陵崗地區(qū)和平原區(qū),平原、低山丘陵和水域面積分別占69%、14%和17%。低山丘陵分布在西南和東北部;平原由北至南為徐淮黃泛平原、里下河平原、濱海平原、沿江平原和太湖平原。主要的土壤類型有水稻土、潮土、砂姜黑土、濱海鹽土、黃棕壤、黃褐土、棕壤等。水田和旱地分別占全省面積的42.46%和23.57%(2006年)。江蘇省是中國重要的糧食生產(chǎn)基地,隨著農(nóng)業(yè)現(xiàn)代化和城市化的迅猛發(fā)展,土壤保護(hù)與利用之間的矛盾比較突出,對土壤空間分布信息的需求非常迫切。

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

    土壤數(shù)據(jù)來源于江蘇省第二次土壤普查資料(1979-1982年)中記錄的典型土壤剖面,主要包括《江蘇土種志》[27]和各市、縣土種志共60本。土壤剖面數(shù)據(jù)主要包括采樣點(diǎn)位置和景觀環(huán)境描述和土壤理化性質(zhì)(主要有SOM、pH值、物理性黏粒和砂粒含量等),篩選出土壤數(shù)據(jù)比較完備的1 519個樣點(diǎn)(圖1)。本研究以表層(0~20 cm)SOM含量為研究對象。對于深度大于20 cm的土層數(shù)據(jù)不做處理,對于小于20 cm的土層,以土層深度為權(quán)重取0~20 cm內(nèi)的所有土層的加權(quán)平均值。SOM含量采用重鉻酸鉀(K2Cr2O7)氧化-滴定法測定。

    年均氣溫和年均降水量數(shù)據(jù),來自中國農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所中國生態(tài)環(huán)境背景層面建造項目完成的柵格數(shù)據(jù)(1 km分辨率),為1980-1999年的逐月平均值計算合成。在ArcGIS軟件中,根據(jù)地形圖、行政區(qū)劃變更信息等資料,將剖面點(diǎn)位置空間化,從氣候柵格數(shù)據(jù)中提取各剖面點(diǎn)的年均氣溫和年均降水?dāng)?shù)據(jù)。

    江蘇省地形以平原為主,研究發(fā)現(xiàn)SOM含量與海拔和坡度等地形因子的相關(guān)性不足0.01,利用逐步回歸建模分析時,無法進(jìn)入模型,可用于SOM空間建模的有效變量較少。對于省級區(qū)域尺度,氣候等條件對土壤屬性的影響較大。鑒于以上背景和原因本文選擇年均氣溫、年均降雨、物理性黏粒含量、土壤pH值等環(huán)境變量進(jìn)行SOM含量的空間建模。

    圖1 土壤樣點(diǎn)分布圖

    1.3 研究方法

    1.3.1 地理加權(quán)回歸

    GWR模型的一般形式如下

    式中(u,v)是第個樣點(diǎn)的空間坐標(biāo);為驗(yàn)證樣點(diǎn)個數(shù);YX是因變量和自變量集X在空間位置(u,v)處的實(shí)測值;0(u,v)為在空間位置(u,v)處的常數(shù)項;系數(shù)(u,v)是連續(xù)函數(shù)(,)在點(diǎn)的值;為符合(0,2)分布的誤差項。通常采用加權(quán)最小二乘法進(jìn)行系數(shù)的局部估計。一般采用高斯(Gaussian)型或雙重平方(bi-square)型空間核函數(shù),利用空間位置(u,v)與周圍一定距離(又稱帶寬)內(nèi)觀測樣點(diǎn)的距離來估算其權(quán)重。本研究采用雙重平方函數(shù)估算權(quán)重,采用黃金分割搜索(golden section search)方法、赤池信息量準(zhǔn)則(akaike information criterion,AIC)[28],確定最優(yōu)帶寬,評價擬合模型的優(yōu)劣。根據(jù)Fotheringham等[11]提出的方法準(zhǔn)則,當(dāng)GWR模型AIC值最小時,帶寬b為最優(yōu)帶寬,此時模型擬合最優(yōu)。詳細(xì)的GWR模型擬合方法參考文獻(xiàn)[11]。如果公式(1)中系數(shù)(,)在空間中任意一點(diǎn)保持不變,則為全局模型,即傳統(tǒng)回歸模型,通常用普通最小二乘法(ordinary least square, OLS)估計參數(shù)。

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

    本研究采用相關(guān)性分析、回歸分析等,研究年均氣溫、年均降雨、物理性黏粒、土壤pH值對江蘇省SOM的影響,采用容差和方差膨脹因子(VIF)對環(huán)境因子與SOM關(guān)系進(jìn)行共線性檢驗(yàn),在SPSS 18.0中完成。在利用GWR進(jìn)行SOM空間建模時,選擇302個樣本作為驗(yàn)證集,1 217個樣本為建模集。為比較不同樣本數(shù)量下GWR建模的精度差異,從建模集中分別隨機(jī)抽取100%(1 217個)、80%(973個)、60%(730個)、40%(486個),20%(243個)的樣點(diǎn)進(jìn)行建模,并選擇最優(yōu)模型進(jìn)行SOM預(yù)測制圖。利用ArcGIS 9.3的地統(tǒng)計模塊中創(chuàng)建子集(create subsets)功能完成不同數(shù)量建模樣本的抽樣。利用GWR4.0軟件進(jìn)行GWR建模與制圖。同時,利用ArcGIS 9.3的地統(tǒng)計模塊以全部建模集(1 217個)進(jìn)行普通克里格插值(ordinary kriging,OK),繪制SOM空間分布圖與GWR預(yù)測制圖結(jié)果做對比。

    在SOM空間預(yù)測制圖時,需要環(huán)境因子(建模的自變量)的空間分布。本研究中物理性黏粒含量和pH數(shù)據(jù)為點(diǎn)狀數(shù)據(jù),利用1 519個樣點(diǎn)進(jìn)行OK插值生成其空間分布,用于SOM空間制圖。綜合考慮氣象數(shù)據(jù)的分辨率、研究區(qū)較大及數(shù)據(jù)運(yùn)算等問題,制圖結(jié)果均采用1 km空間分辨率。

    1.3.3 模型評價

    選取平均誤差(mean error,ME)、平均絕對誤差(mean absolute error,MAE)、均方根誤差(root mean of squared error,RMSE)和決定系數(shù)(determination coefficient,2)等指標(biāo)評價建模結(jié)果,其計算公式如下

    2 結(jié)果與分析

    2.1 SOM主要影響因素分析

    江蘇省SOM質(zhì)量分?jǐn)?shù)范圍較大(= 1 519),為1.3~52.4 g/kg,極差為51.1 g/kg;平均為16.55 g/kg,標(biāo)準(zhǔn)差為8.49 g/kg(表1)。變異系數(shù)為51.36%,屬于中等變異強(qiáng)度。偏度和峰度系數(shù)分別為1.12和1.75,經(jīng)自然對數(shù)變換后符合正態(tài)分布。不同樣點(diǎn)數(shù)量的建模集的SOM頻率分布均符合對數(shù)正態(tài)分布。

    相關(guān)分析表明,江蘇省SOM與年均氣溫極顯著正相關(guān),相關(guān)系數(shù)為0.47(<0.001,=1 519),以年均降雨量為控制變量時二者的偏相關(guān)系數(shù)為0.21(<0.001,= 1 519),即在溫度較高的長江以南地區(qū)SOM含量較高,溫度較低的蘇北地區(qū)含量較低。這與周濤等[9]、許信旺等[29]研究結(jié)果一致。SOM與年均降雨量極顯著正相關(guān),相關(guān)系數(shù)為0.44(<0.001,=1 519),以年均氣溫為控制變量時二者的相關(guān)性降低,偏相關(guān)系數(shù)為0.11(<0.001,=1 519)。這表明SOM與降雨量的相關(guān)性是由年均降雨與年均氣溫的相關(guān)性(=0.82,<0.001)引起的。SOM與土壤物理性黏粒和砂粒含量相關(guān)系數(shù)為0.46和?0.45(<0.001,=1 519)。SOM與pH值的相關(guān)系數(shù)為?0.30(<0.001,=1 519)。線性回歸分析表明年均氣溫、年均降雨、土壤質(zhì)地和pH值對SOM含量(=1 519)變異的綜合解釋能力為45.3%(0.001)?;貧w模型中的各因子的容差介于0.31~0.98,均大于0.1,VIF介于0.02~3.25,均不大于7.5,表明自變量之間不存在高度共線性[30]。

    表1 江蘇省SOM含量統(tǒng)計值

    注:括號內(nèi)為自然對數(shù)變換值。

    Note: values in the brackets were logarithm transformed.

    2.2 SOM空間自相關(guān)分析

    利用ArcGIS軟件空間統(tǒng)計功能分析不同建模樣本集的SOM的全局和局部空間自相關(guān)性。在ArcGIS中計算全局Moran’s I值時,在0~100 km范圍內(nèi)按照10 km間隔設(shè)定不同的空間距離,計算全局Moran’s I值。圖2為不同建模集的全局Moran’s I值,結(jié)果表明不同建模樣本集的SOM在不同的空間尺度上均存在不同程度的極顯著空間自相關(guān)性,全局Moran’s I值介于0.25~0.61(<0.001)。隨著空間距離的增大,不同建模集的SOM的全局Moran’s I值均逐漸減小,表明其空間自相關(guān)性隨著空間距離的增大而遞減。在相同的空間距離上,除建模集為=486時的全局Moran’s I值稍小外,其他建模集的SOM的全局Moran’s I值總體上相差不大。

    圖3為建模集為=1 217時的局部空間自相關(guān)聚類圖。由圖知,江蘇省SOM含量存在空間聚集特征,“高-高”(high-high)聚集區(qū)(SOM含量較高且空間差異較小的樣點(diǎn)聚集分布)主要分布在蘇中和蘇南地區(qū),“低-低”(low-low)聚集區(qū)(SOM含量較低且空間差異較小的樣點(diǎn)聚集分布)主要分布在蘇北地區(qū);“高-低”(high-low)聚集區(qū)(SOM含量較高的樣點(diǎn)被含量較低的樣點(diǎn)包圍,且含量的空間差異顯著)主要分布在蘇北地區(qū),“低-高”(low-high)聚集區(qū)(SOM含量較低的樣點(diǎn)被含量較高的樣點(diǎn)包圍,且含量的空間差異顯著)主要分布在蘇南地區(qū)。不同建模集的局部空間自相關(guān)聚類圖結(jié)果相似,沒有一一作圖展示。以上分析表明不同建模集的SOM在不同的空間尺度均存在顯著的全局和局部空間自相關(guān)性,可以采用GWR模型進(jìn)行研究區(qū)SOM建模[17]。

    注:顯著性為P<0.001。

    注:圖例代表不同聚類聚集區(qū),下同。

    2.3 基于GWR的SOM建模分析

    將不同建模集的SOM經(jīng)過自然對數(shù)變換后,分別進(jìn)行GWR建模和全局建模。為了比較不同自變量在建模中的貢獻(xiàn)大小,首先在建模前將各自變量進(jìn)行標(biāo)準(zhǔn)化變換到同一量綱,然后通過其系數(shù)來對比各自變量在SOM建模中的作用。

    表2為不同建模集的SOM的GWR建模和全局建模結(jié)果。一般認(rèn)為GWR模型和全局模型的AIC相差大于3,即使考慮到GWR模型的復(fù)雜性,其也明顯優(yōu)于全局模型。

    表2 不同樣點(diǎn)數(shù)量的GWR建模和全局建模結(jié)果比較

    表2顯示對于不同建模集的GWR模型的AIC均比相應(yīng)的全局模型要小,相差56.08~360.19之間,表明GWR模型明顯優(yōu)于全局模型。不同建模集的GWR模型的調(diào)整決定系數(shù)(2adj)均大于相應(yīng)的全局模型,其2adj較相應(yīng)的全局模型增加0.15~0.20之間,表明對于相同的建模集,GWR模型對SOM變異的解釋能力大于全局模型。對于不同建模集SOM的GWR模型,2adj介于0.54~0.62,總體上隨著建模樣本數(shù)量的減少GWR模型對SOM的解釋能力稍有差異,但差異不大。這說明對于江蘇省SOM的建模估計,采用較少樣點(diǎn)數(shù)進(jìn)行GWR建模時,也可滿足一定的精度需求(如建模樣本= 486和= 1 217的模型擬合度較接近)。對于不同建模集SOM的全局模型,其2adj介于0.39~0.44之間。

    不同建模集的GWR模型的殘差平方和(residual sum of squares,RSS)均比相應(yīng)的全局模型減少17.40~76.67之間。GWR模型和全局模型殘差的方差分析結(jié)果(表3)顯示,不同建模集的GWR模型的殘差較全局模型的殘差顯著降低,其檢驗(yàn)值均達(dá)到0.01極顯著性水平,表明GWR模型的擬合效果較全局模型有顯著改善。上述這些結(jié)果表明,不同建模集的GWR建模結(jié)果均優(yōu)于全局建模,且當(dāng)建模樣點(diǎn)數(shù)較少時GWR模型也能有較好的建模效果。

    表4為不同樣點(diǎn)數(shù)量的GWR模型回歸系數(shù)和全局模型回歸系數(shù)。結(jié)果顯示,在不同樣點(diǎn)數(shù)量的全局建模中(除樣點(diǎn)數(shù)量= 243),同一環(huán)境因子對SOM預(yù)測的作用相似,如各全局模型中黏粒的回歸系數(shù)介于0.24~0.25,相差較??;不同變量對SOM的作用大小依次為:黏粒> 年均溫> 年均降雨> 土壤pH值。在不同樣點(diǎn)數(shù)量的GWR建模中(除樣點(diǎn)數(shù)量= 243),從回歸系數(shù)的中值和平均值來看,各環(huán)境因子對SOM的作用大小規(guī)律與全局建模一致。當(dāng)建模樣點(diǎn)數(shù)量較少(= 243時),環(huán)境因子與SOM之間的關(guān)系不能完全地表達(dá)和刻畫。GWR模型中,各環(huán)境因子的回歸系數(shù)存在不同幅度的變化范圍,同一環(huán)境因子的回歸系數(shù)變化范圍相差不大,其中黏粒的回歸系數(shù)變化范圍最穩(wěn)定且系數(shù)均為正(表4),體現(xiàn)了黏粒含量對SOM含量的正相關(guān)。

    表3 GWR模型和全局建模殘差的方差分析結(jié)果

    注:**為顯著性< 0.01,下同。

    Note: **is significant difference at 0.01 level, the same below.

    圖4為GWR模型回歸系數(shù)(建模樣點(diǎn)= 1 217),各環(huán)境因子對SOM的作用存在顯著的空間差異,即SOM與這些環(huán)境因子之間的關(guān)系存在空間非平穩(wěn)性[19]。各因子對SOM的影響程度可通過對應(yīng)的回歸系數(shù)來解釋,回歸系數(shù)的絕對值大小反映該因子對SOM影響的強(qiáng)烈程度,正反映因子對SOM的正影響或正相關(guān),反之亦然[13, 30]。如,黏粒的回歸系數(shù)變化范圍比其他因素回歸系數(shù)的變化范圍較小,表明在研究區(qū)內(nèi)黏粒對SOM的影響相比其他因素較穩(wěn)定。黏粒對SOM的影響在江蘇省東北部和西南部最小,在中南部最大,且均表現(xiàn)為正的影響(圖4a),這與趙明松等[10]在該地區(qū)的研究結(jié)果相似。年均溫對SOM的作用在大部分地區(qū)表現(xiàn)為正相關(guān)(圖4b);年均降雨對SOM的作用在江蘇省東部地區(qū)表現(xiàn)為負(fù)相關(guān),在西部表現(xiàn)為正相關(guān)。雖然研究區(qū)內(nèi)年均溫和年均降雨顯著正相關(guān),但是對于SOM的作用強(qiáng)度在空間上顯著差異,相應(yīng)系數(shù)的空間格局差異較大。土壤pH值對SOM的作用,在西北-東南方向上向南北兩側(cè)逐漸變化,總體上由負(fù)相關(guān)轉(zhuǎn)變?yōu)檎嚓P(guān)。

    表4 不同樣點(diǎn)數(shù)量的GWR模型回歸系數(shù)

    注:*顯著性< 0.05,下同。

    Note: *, ** Siginificant difference at 0.05 level, the same below.

    圖4 GWR建?;貧w系數(shù)空間分布圖(n=1 217)

    2.4 基于GWR的SOM預(yù)測制圖及模型評價

    圖5為江蘇省SOM含量空間分布預(yù)測圖,其中a~e為不同樣點(diǎn)數(shù)量的GWR建模預(yù)測圖,f為采用全部建模樣點(diǎn)(=1 217)的OK插值圖。

    圖5 SOM含量空間分布預(yù)測圖

    SOM的OK插值的變異函數(shù)采用球狀模型擬合,塊金值為0.13,偏基臺值為0.10,變程為75 km,步長為10 km。結(jié)果顯示,不同樣點(diǎn)數(shù)量的GWR預(yù)測結(jié)果的范圍介于3.38~45.08 g/kg,平均值介于14.94~15.34 g/kg,標(biāo)準(zhǔn)差介于5.99~7.15 g/kg,均比建模樣點(diǎn)的相應(yīng)統(tǒng)計值稍小。不同樣點(diǎn)數(shù)量的GWR預(yù)測結(jié)果的空間分布格局基本一致:SOM含量的最高值主要分布在太湖平原和里下河平原,從蘇南到蘇北地區(qū)逐漸降低;在研究區(qū)西南部的部分地區(qū)SOM的空間分布細(xì)節(jié)存在細(xì)微差異(圖5)。總體上,不同樣點(diǎn)數(shù)量的GWR預(yù)測結(jié)果與采用全部建模樣點(diǎn)的OK預(yù)測結(jié)果具有較相似的空間分布格局。

    表5為SOM建模的獨(dú)立驗(yàn)證結(jié)果。不同樣點(diǎn)數(shù)量(除= 243時)的GWR建模的驗(yàn)證結(jié)果差異不大:RMSE介于5.56~5.75 g/kg,MAE介于3.87~4.05 g/kg,均小于全部建模樣點(diǎn)的OK驗(yàn)證結(jié)果;2介于0.48~0.52,均高于OK插值的驗(yàn)證結(jié)果。不同樣點(diǎn)數(shù)量的GWR驗(yàn)證結(jié)果的ME介于?0.58~?0.22 g/kg,表明GWR建模結(jié)果總體上存在不同程度的高估現(xiàn)象。樣點(diǎn)數(shù)=243時GWR建模的驗(yàn)證結(jié)果(RMSE=6.31,ME=?0.22,MAE=4.47)較其他稍差,但其預(yù)測的空間分布格局與其他結(jié)果較相似。上述結(jié)果表明,在江蘇省利用GWR進(jìn)行SOM空間預(yù)測制圖時,樣點(diǎn)數(shù)量對建模精度影響不大,較少的樣點(diǎn)數(shù)量也可以獲得較好的制圖結(jié)果。與OK插值方法比較發(fā)現(xiàn),在樣點(diǎn)數(shù)最大時,GWR建模精度提高幅度有限,2僅提高了0.06。但從另一個方面看,如果樣點(diǎn)較少時(如486個)采用GWR建模,能夠獲得和較多樣點(diǎn)的OK插值結(jié)果的相似精度,對減少制圖成本也有一定的貢獻(xiàn)。

    表5 不同樣本數(shù)量的GWR建模的驗(yàn)證結(jié)果

    圖6為不同樣本數(shù)量的GWR模型殘差的全局Moran’I指數(shù)。GWR模型殘差的全局Moran’s I值介于0.01~0.12(>0.1),隨著空間距離的增大其值均逐漸減小,表明在不同的空間尺度上不同樣點(diǎn)數(shù)量的GWR模型殘差均不存在顯著空間自相關(guān)性。在相同的空間距離上,GWR殘差的全局Moran’s I值總體上相差不大。圖7為GWR模型殘差(= 1 217)的局部空間自相關(guān)聚類圖。由圖知,絕大部分樣點(diǎn)的建模均不存在顯著局部空間自相關(guān)性,僅有少量樣點(diǎn)(= 45)的殘差存在著“高-高”(殘差較高且空間差異較小的樣點(diǎn)聚集分布區(qū))、“低-低”和“高-低”聚類。不同建模集的GWR殘差的局部空間自相關(guān)聚類圖結(jié)果相似,沒有一一作圖展示。不同樣點(diǎn)數(shù)量GWR建模的殘差空間自相關(guān)性較弱,也間接表明GWR擬合的模型較優(yōu)。

    注:不顯著(P=0.1)。

    圖7 GWR模型殘差局部空間聚類圖(n=1 217)

    3 結(jié) 論

    1)江蘇省SOM含量在不同空間尺度上存在極顯著的全局和局部空間自相關(guān)性。SOM含量空間分布以空間聚集特征為主,“高-高”聚集區(qū)主要分布在蘇中和蘇南地區(qū),“低-低”聚集區(qū)主要分布在蘇北地區(qū)。不同建模集的全局自相關(guān)性和局部空間自相關(guān)聚類圖結(jié)果相似。

    2)GWR建模結(jié)果均優(yōu)于全局建模,對SOM空間變異的解釋能力更高,不同建模集的GWR的2adj較全局建模均有明顯提高,其AIC和RSS均比全局模型有大幅降低,擬合效果有較大提升。隨著建模樣本數(shù)量的減少GWR模型對SOM的解釋能力差異不大。GWR建模的殘差在不同的空間尺度上均不存在空間自相關(guān)性。

    3)利用GWR模型進(jìn)行SOM空間預(yù)測制圖時,建模樣點(diǎn)數(shù)量對建模精度影響不大,不同樣點(diǎn)數(shù)量的GWR建模結(jié)果與采用全部樣點(diǎn)的OK插值結(jié)果具有相似的空間格局。樣點(diǎn)數(shù)較少時,在大區(qū)域利用GWR模型進(jìn)行SOM預(yù)測制圖,也可以達(dá)到一定的精度要求。

    數(shù)字土壤制圖過程中,一般要求環(huán)境因子為空間柵格形式。本文中物理性黏粒含量和pH值為點(diǎn)狀數(shù)據(jù),在預(yù)測制圖中首先通過插值其空間分布,再用于SOM預(yù)測制圖,可能增加二次誤差,對最終的制圖結(jié)果可能會有一定程度的影響,這些有待進(jìn)一步的深入研究。

    [1] Yadav V, Malanson G. Progress in soil organic matter research: Litter decomposition, modelling, monitoring and sequestration [J]. Progress in Physical Geography, 2007, 31(2): 131-154.

    [2] Davidson E A, Trumbore S E, Amundson R. Biogeochemistry: Soil warming and organic carbon content[J]. Nature, 2000, 408(6814): 789-790.

    [3] McBratney A B, Santos M L M, Minasny B. On digital soil mapping [J]. Geoderma, 2003, 117: 3-52.

    [4] Thompson J A, Pena-Yewtukhiw E M, Grove J H. Soil-landscape modeling across a physiographic region: Topographic patterns and model transportability [J]. Geoderma, 2006, 133: 57-70.

    [5] Qi F, Zhu A X, Harrower M, et al. Fuzzy soil mapping based on prototype category theory[J]. Geoderma, 2006, 136: 774-787.

    [6] Zhu A X, Qi F, Moore A, et al. Prediction of soil properties using fuzzy membership values [J]. Geoderma, 2010, 158: 199-206.

    [7] de Bruin S, Stein A. Soil-landscape modelling using fuzzy c-means clustering of attribute data derived from a digital elevation model (DEM) [J]. Geoderma, 1998, 83: 17-33.

    [8] 朱阿興,楊琳,樊乃卿,等. 數(shù)字土壤制圖研究綜述與展望[J]. 地理科學(xué)進(jìn)展,2018,37(1):66-78.

    Zhu Axing, Yang Lin, Fan Naiqing, et al. The review and outlook of digital soil mapping[J]. Progress in Geography, 2018, 37(1): 66-78. (in Chinese with English abstract)

    [9] 周濤,史培軍,王紹強(qiáng). 氣候變化及人類活動對中國土壤有機(jī)碳儲量的影響[J]. 地理學(xué)報,2003,58(5):727-734.

    Zhou Tao, Shi Peijun, Wang Shaoqiang. Impacts of climate and human activities on soil carbon storage in China[J]. Acta Geographica Sinica, 2003, 58(5): 727-734. (in Chinese with English abstract)

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

    Zhao Mingsong, Zhang Ganlin, Li Decheng, et al. Variability of soil organic matter and its main factors in Jiangsu Province[J]. Acta Ecologica Sinica, 2013, 33(16): 5058-5066. (in Chinese with English abstract)

    [11] Fotheringham A S, Brunsdon C, Charlton M. Geographically Weighted Regression—the Analysis of Spatially Varying Relationships[M]. Chichester.UK: John Wiley & Sons Inc., 2002: 42-46.

    [12] 王遠(yuǎn)飛,何洪林. 空間數(shù)據(jù)分析方法[M]. 北京:科學(xué)出版社,2007:132-141.

    [13] 瞿明凱,李衛(wèi)東,張傳榮,等. 地理加權(quán)回歸及其在土壤和環(huán)境科學(xué)上的應(yīng)用前景[J]. 土壤,2014,46(1):15-22.

    Qu Mingkai, Li Weidong, Zhang Chuanrong, et al. Geographically weighted regression and its application prospect in soil and environmental sciences[J]. Soils, 2014, 46(1): 15-22. (in Chinese with English abstract)

    [14] 呂萍,甄輝. 基于GWR模型的北京市住宅用地價格影響因素及其空間規(guī)律研究[J]. 經(jīng)濟(jì)地理,2010,30(3):472-478.

    Lü Ping, Zhen Hui. Affecting factors research of Beijing residential land price based on GWR model[J]. Economic Geography, 2010, 30(3): 472-478. (in Chinese with English abstract)

    [15] 湯慶園,徐偉,艾福利. 基于地理加權(quán)回歸的上海市房價空間分異及其影響因子研究[J]. 經(jīng)濟(jì)地理,2012,32(2):52-58.

    Tang Qingyuan, Xu Wei, Ai Fuli. A GWR based study on spatial pattern and structural determinants of Shanghai’s housing price[J]. Economic Geography, 2012, 32(2): 52-58. (in Chinese with English abstract)

    [16] 劉衛(wèi)東,劉紅光,范曉梅,等. 地區(qū)間貿(mào)易流量的產(chǎn)業(yè):空間模型構(gòu)建與應(yīng)用[J]. 地理學(xué)報,2012,67(2):147-156.

    Liu Weidong, Liu Hongguang, Fan Xiaomei, et al. Sector-specific spatial statistic model for estimating inter-regional Trade tflows: A case study of agricultural, chemical and electronic sectors in China[J]. Acta Geographica Sinica, 2012, 67(2): 147-156. (in Chinese with English abstract)

    [17] 劉瓊峰,李明德,段建南,等. 農(nóng)田土壤鉛、鎘含量影響因素地理加權(quán)回歸模型分析[J]. 農(nóng)業(yè)工程學(xué)報,2013,29(3):225-234.

    Liu Qiongfeng, Li Mingde, Duan Jiannan, et al. Analysis on influence factors of soil Pb and Cd in agricultural soil of Changsha suburb based on geographically weighted regression model [J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2013, 29(3): 225-234. (in Chinese with English abstract)

    [18] 霍霄妮,李紅,孫丹峰,等. 北京耕地土壤重金屬空間自回歸模型及影響因素[J]. 農(nóng)業(yè)工程學(xué)報,2010,26(5):78-82.

    Huo Xiaoni, Li Hong, Sun Danfeng, et al. Spatial autogression model for heavy metals in cultivated soils of Beijing[J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2010, 26(5): 78-82. (in Chinese with English abstract)

    [19] 李錦芬,瞿明凱,黃標(biāo),等. 區(qū)域土壤CEC與相關(guān)控制因子的空間非平穩(wěn)關(guān)系評估[J]. 土壤學(xué)報,2017,54(3):639-647.

    Li Jinfen, Qu Mingkai, Huang Biao, et al. Spatially non-stationary relationships between cation exchange capacity and related control factors[J]. Acta Pedologica Sinica, 2017, 54(3): 639-647. (in Chinese with English abstract)

    [20] 王景雷,康紹忠,孫景生,等. 基于PCA和GWR的作物需水量空間分布估算[J]. 科學(xué)通報,2013,58(12):1131-1139.

    Wang Jinglei, Kang Shaozhong, Sun Jingsheng, et al. Estimation of crop water requirement based on principal component analysis and geographically weighted regression[J]. China Science Bulletin, 2013, 58, 1131-1139. (in Chinese with English abstract)

    [21] Song X D, Brus D J, Liu F, et al. Mapping soil organic carbon content by geographically weighted regression: A case study in the Heihe River Basin, China[J]. Geoderma, 2016, 261: 11-22.

    [22] Mishra U, Lal R, Liu D S, et al. Predicting the spatial variation of the soil organic carbon pool at a regional scale[J]. Soil Science Society of America Journal, 2010, 74(3): 906-914.

    [23] Yang S H, Liu F, Song X D, et al. Mapping topsoil electrical conductivity by a mixed geographically weighted regression kriging: A case study in the Heihe River Basin, northwest China[J]. Ecological Indicators, 2019, 102: 252-264.

    [24] Zhang C S, Tang Y, Xu X L, et al. Towards spatial geochemical modelling: Use of geographically weighted regression for mapping soil organic carbon contents in Ireland[J]. Applied Geochemistry, 2011, 26(7): 1239-1248.

    [25] Zeng C, Yang L, Zhu A X, et al. Mapping soil organic matter concentration at different scales using a mixed geographically weighted regression method[J]. Geoderma, 2016, 281: 69-82.

    [26] 郭龍,張海濤,陳家贏,等. 基于協(xié)同克里格插值和地理加權(quán)回歸模型的土壤屬性空間預(yù)測比較[J]. 土壤學(xué)報,2012,49(5):1037-1042.

    Guo Long, Zhang Haitao, Chen Jiaying, et al. Comparison between cokring model and geographically weighted regression model in spatial prediction of soil attributes[J]. Acta Pedologica Sinica, 2012, 49(5): 1037-1042. (in Chinese with English abstract)

    [27] 江蘇省土壤普查辦公室. 江蘇土種志[M]. 南京:江蘇科學(xué)技術(shù)出版社,1996:5-105.

    [28] Akaike H. A new look at the statistical model identification[J]. IEEE Transactions on Automatic Control, 1974, 19(6): 716-723.

    [29] 許信旺,潘根興,汪艷林,等. 中國農(nóng)田耕層土壤有機(jī)碳變化特征及控制因素[J]. 地理研究,2009,28(3):601-612.

    Xu Xinwang, Pan Genxing, Wang Yanlin, et al. Research of changing characteristics and control factors of farmland topsoil organic carbon in China[J]. Geographical Research, 2009, 28(3): 601-612. (in Chinese with English abstract)

    [30] 許爾琪. 基于地理加權(quán)回歸的石漠化影響因子分布研究[J]. 資源科學(xué),2017,39(10):1975-1988.

    Xu Erqi. Spatial variation in drivers of karst rocky desertification based on geographically weighted regression model[J]. Resources Science, 2017, 39(10): 1975-1988. (in Chinese with English abstract)

    Spatial modeling of soil organic matter over low relief areas based on geographically weighted regression

    Zhao Mingsong1,2, Liu Binyin1, Lu Hongliang1, Li Decheng2, Zhang Ganlin2※

    (1.232001; 2.210008)

    Accurate estimates of the spatial variability of soil organic matter (SOM) are necessary to properly evaluate climatic chagne, soil carbon sequestration potential and soil fertility. In plains and gently undulating terrains, soil spatial variability is not closely related to relief, and thus digital soil mapping (DSM) methods based on soil-landscape relationships often fail in these areas. Therefore, different predictors or methods are needed for DSM in plains. In provincial regional scale, climatic factors influence spatial distribution of soil properties. For this research, Jiangsu Province was selected as example and mean annual temperature (MAT), mean annual precipitation (MAP), physical clay content, and soil pH were selected for SOM spatial modeling using geographically weighted regression (GWR). The SOM content in the surface layer (0-20cm) of 1 519 typical soil profiles of the Second National Soil Survey in Jiangsu Province were collected. 1 217 samples were selected as the modeling set and 302 were the validation set. Fristly, 100% (1 217), 80% (973), 60% (730), 40% (486), and 20% (243) samples were randomly selected from the modeling set, and global and local spatial autocorrelation of SOM content were analyzed at different spatial scales using spatial statistics tools in ArcGIS. Secondly, comparison of the accuracy between GWR model and the global regression model under the different sampling size was conducted. Akaike information criterion (AIC), residual sum of squares (RSS) and adjustment determination coefficient (2adj) were used modeling comparison. Thirdly, the optimal model was selected for mapping SOM spatial prediction. Independent validation was used for model evaluation, using four indices: mean error (ME), mean absolute error (MAE) and root mean of squared error (RMSE), and determination coefficient (2). Results show that: 1) There was a significant spatial autocorrelation of SOM content in Jiangsu Province at different spatial scales. The clustering pattern of global and local spatial autocorrelation of modeling set with different sampling size were similar. The global Moran’s I ranged from 0.25 to 0.61 (<0.001). The spatial distribution of SOM content was mainly characterized by spatial clustering pattern. The “high-high” clustering areas were mainly distributed in the central and south of Jiangsu, and the “l(fā)ow-low” clustering areas were mainly distributed in the north of Jiangsu. 2) The modeling results of GWR were better than the global regression modeling, and the residuals had no spatial autocorrelation at different spatial scales. The2adjof GWR in different modeling sets was increased by 0.15 to 0.20 compared with the global model. The AIC and RSS were significantly lower than the global model, which were decreased by 56.08 to 360.19 and 17.40 to 76.67 respectively. There were slight difference between GWR models with different sampling size. 3) The number of modeling samples (except for the number of modeling samples was 243) had little effect on the accuracy of prediction and mapping results of GWR, the RMSE was between 5.56 and 5.75 g/kg, MAE was between 3.87 and 4.05 g/kg and2was between 0.48 and 0.52. The results were all better than the validation result of Ordinary Kriging using all modeling sampling points. This study can provide reference for SOM modeling and mapping in large and low relief areas with sparse samples.

    soils; organic matter; models; geographically weighted regression; digital soil mapping; low relief areas; Jiangsu Province

    趙明松,劉斌寅,盧宏亮,李德成,張甘霖. 基于地理加權(quán)回歸的地形平緩區(qū)土壤有機(jī)質(zhì)空間建模[J]. 農(nóng)業(yè)工程學(xué)報,2019,35(20):102-110.doi:10.11975/j.issn.1002-6819.2019.20.013 http://www.tcsae.org

    Zhao Mingsong, Liu Binyin, Lu Hongliang, Li Decheng, Zhang Ganlin. Spatial modeling of soil organic matter over low relief areas based on geographically weighted regression[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(20): 102-110. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.2019.20.013 http://www.tcsae.org

    2019-03-26

    2019-06-23

    國家自然科學(xué)基金(41501226);土壤與農(nóng)業(yè)可持續(xù)發(fā)展國家重點(diǎn)實(shí)驗(yàn)室開發(fā)基金(Y412201431);安徽省高校自然科學(xué)研究項目(KJ2015A034)

    趙明松,博士,副教授,主要從事數(shù)字土壤制圖和空間變異研究。Email:zhaomingsonggis@163.com

    張甘霖,博士,研究員,博士生導(dǎo)師,主要從事土壤發(fā)生分類、土壤地理研究。Email:glzhang@issas.ac.cn

    10.11975/j.issn.1002-6819.2019.20.013

    S159.9

    A

    1002-6819(2019)-20-0102-09

    猜你喜歡
    樣點(diǎn)回歸系數(shù)制圖
    小麥條銹病田間為害損失的初步分析
    湖北植保(2022年4期)2022-08-23 10:51:52
    基于空間模擬退火算法的最優(yōu)土壤采樣尺度選擇研究①
    土壤(2021年1期)2021-03-23 07:29:06
    無聲手槍如何消音?
    多元線性回歸的估值漂移及其判定方法
    基于分融策略的土壤采樣設(shè)計方法*
    電導(dǎo)法協(xié)同Logistic方程進(jìn)行6種蘋果砧木抗寒性的比較
    二向反射模型在土地覆被制圖中的應(yīng)用
    多元線性模型中回歸系數(shù)矩陣的可估函數(shù)和協(xié)方差陣的同時Bayes估計及優(yōu)良性
    工程制圖課程教學(xué)改革探析
    建筑工程制圖與識圖專業(yè)人才培養(yǎng)的探討
    河南科技(2014年3期)2014-02-27 14:06:14
    色婷婷久久久亚洲欧美| 国产男女超爽视频在线观看| 日韩精品有码人妻一区| av天堂久久9| 亚洲国产精品成人久久小说| 久久久精品免费免费高清| 欧美xxⅹ黑人| 国产精品三级大全| 在线免费观看不下载黄p国产| 久久精品熟女亚洲av麻豆精品| 亚洲成人av在线免费| 亚洲欧洲国产日韩| 亚洲高清免费不卡视频| 午夜福利视频在线观看免费| 亚洲欧美精品自产自拍| 有码 亚洲区| 中文乱码字字幕精品一区二区三区| tube8黄色片| 99国产综合亚洲精品| 精品国产一区二区久久| 亚洲精品国产av成人精品| 亚洲图色成人| 日本-黄色视频高清免费观看| 七月丁香在线播放| 日韩免费高清中文字幕av| 精品亚洲乱码少妇综合久久| videosex国产| 99热网站在线观看| 欧美国产精品va在线观看不卡| 精品人妻一区二区三区麻豆| 精品国产一区二区三区久久久樱花| 中文字幕人妻熟女乱码| 9191精品国产免费久久| 又黄又爽又刺激的免费视频.| 美女脱内裤让男人舔精品视频| 日本wwww免费看| 少妇猛男粗大的猛烈进出视频| 亚洲精品中文字幕在线视频| 久久久国产精品麻豆| 欧美日韩精品成人综合77777| 免费大片18禁| 成人二区视频| 九草在线视频观看| 我的女老师完整版在线观看| 亚洲精品国产av成人精品| av在线app专区| 亚洲成人av在线免费| 国产麻豆69| 卡戴珊不雅视频在线播放| 久久久久久久国产电影| 你懂的网址亚洲精品在线观看| 桃花免费在线播放| 我要看黄色一级片免费的| 日本91视频免费播放| 国产精品久久久久久久电影| 中国三级夫妇交换| 国产精品嫩草影院av在线观看| 国产色爽女视频免费观看| 香蕉丝袜av| 日韩三级伦理在线观看| 97人妻天天添夜夜摸| 人妻 亚洲 视频| 一级,二级,三级黄色视频| www.熟女人妻精品国产 | 成人手机av| 天堂8中文在线网| 999精品在线视频| 亚洲久久久国产精品| 波多野结衣一区麻豆| 91在线精品国自产拍蜜月| 亚洲国产成人一精品久久久| av有码第一页| 桃花免费在线播放| 亚洲人与动物交配视频| 国产成人一区二区在线| 国产麻豆69| 一级毛片我不卡| 欧美精品高潮呻吟av久久| 国产免费福利视频在线观看| 黑人欧美特级aaaaaa片| 亚洲图色成人| 亚洲精品aⅴ在线观看| 大片电影免费在线观看免费| 国产乱来视频区| 免费在线观看完整版高清| 免费人妻精品一区二区三区视频| 美女内射精品一级片tv| 亚洲精品久久成人aⅴ小说| 亚洲激情五月婷婷啪啪| 亚洲精品久久午夜乱码| 精品99又大又爽又粗少妇毛片| 免费av不卡在线播放| 最近最新中文字幕免费大全7| 欧美人与性动交α欧美软件 | 丰满迷人的少妇在线观看| 成人二区视频| 精品国产一区二区三区久久久樱花| 各种免费的搞黄视频| 五月天丁香电影| 欧美变态另类bdsm刘玥| 日韩中字成人| 蜜臀久久99精品久久宅男| 午夜免费鲁丝| 国产高清三级在线| 男女啪啪激烈高潮av片| 一区二区三区四区激情视频| 91成人精品电影| 亚洲成人av在线免费| 色吧在线观看| 熟女人妻精品中文字幕| 看免费成人av毛片| 最近最新中文字幕免费大全7| 在线 av 中文字幕| 亚洲精品美女久久久久99蜜臀 | 久久久久国产精品人妻一区二区| 亚洲成av片中文字幕在线观看 | 一边亲一边摸免费视频| 青春草国产在线视频| 国产一区二区三区综合在线观看 | 久久这里有精品视频免费| 2022亚洲国产成人精品| 9191精品国产免费久久| 亚洲三级黄色毛片| 婷婷色综合大香蕉| 亚洲精品日本国产第一区| 久久久欧美国产精品| 国产 一区精品| 大片免费播放器 马上看| 国产成人精品一,二区| 日韩人妻精品一区2区三区| 美女xxoo啪啪120秒动态图| 飞空精品影院首页| 天天操日日干夜夜撸| 亚洲国产精品一区三区| videos熟女内射| 精品福利永久在线观看| 欧美精品av麻豆av| 少妇人妻精品综合一区二区| 丝袜脚勾引网站| 婷婷色综合大香蕉| 国产精品秋霞免费鲁丝片| 一本—道久久a久久精品蜜桃钙片| 美女xxoo啪啪120秒动态图| 一本久久精品| 如何舔出高潮| 久久亚洲国产成人精品v| 人人妻人人澡人人爽人人夜夜| 在线天堂最新版资源| 亚洲av国产av综合av卡| 丝袜人妻中文字幕| 国产一区二区在线观看av| 久久ye,这里只有精品| 一区二区三区乱码不卡18| 51国产日韩欧美| 国产深夜福利视频在线观看| 久久精品国产鲁丝片午夜精品| 另类亚洲欧美激情| 女性生殖器流出的白浆| 黄片无遮挡物在线观看| 色婷婷久久久亚洲欧美| 亚洲精品一二三| 少妇人妻 视频| 国产老妇伦熟女老妇高清| 日韩成人伦理影院| 美女国产视频在线观看| 国产精品国产三级国产专区5o| 亚洲在久久综合| www.色视频.com| 最近2019中文字幕mv第一页| 精品国产一区二区三区久久久樱花| 国产精品三级大全| 九九爱精品视频在线观看| 成人毛片60女人毛片免费| 视频中文字幕在线观看| 精品国产一区二区三区四区第35| av播播在线观看一区| 嫩草影院入口| 国产毛片在线视频| 少妇的逼好多水| 亚洲av综合色区一区| 熟女电影av网| 亚洲 欧美一区二区三区| 欧美精品高潮呻吟av久久| 蜜臀久久99精品久久宅男| 18禁在线无遮挡免费观看视频| h视频一区二区三区| 免费高清在线观看视频在线观看| 夜夜骑夜夜射夜夜干| 另类亚洲欧美激情| 五月玫瑰六月丁香| 自线自在国产av| 国产永久视频网站| 成人无遮挡网站| 国产又色又爽无遮挡免| 日本vs欧美在线观看视频| 欧美xxxx性猛交bbbb| 成人免费观看视频高清| 久久综合国产亚洲精品| 日本欧美视频一区| 少妇人妻 视频| 51国产日韩欧美| 国产成人精品在线电影| 欧美日韩视频高清一区二区三区二| 老女人水多毛片| 大香蕉久久成人网| 又黄又粗又硬又大视频| 欧美人与善性xxx| 各种免费的搞黄视频| 三级国产精品片| 99视频精品全部免费 在线| 中文字幕最新亚洲高清| 成人综合一区亚洲| 日韩伦理黄色片| a 毛片基地| 成人亚洲精品一区在线观看| 男女国产视频网站| 精品一区二区三卡| 99久久中文字幕三级久久日本| 女性被躁到高潮视频| 成年女人在线观看亚洲视频| 欧美日韩一区二区视频在线观看视频在线| 亚洲熟女精品中文字幕| 观看美女的网站| 亚洲美女搞黄在线观看| 久久精品人人爽人人爽视色| 亚洲成国产人片在线观看| 狂野欧美激情性bbbbbb| 日本欧美视频一区| 亚洲第一区二区三区不卡| 综合色丁香网| 午夜免费观看性视频| 日本wwww免费看| 国产日韩欧美视频二区| 成人毛片60女人毛片免费| 日日摸夜夜添夜夜爱| 两个人免费观看高清视频| 国产欧美日韩综合在线一区二区| 欧美丝袜亚洲另类| 国产黄色免费在线视频| 大香蕉97超碰在线| 最近中文字幕2019免费版| 免费大片18禁| 久久精品aⅴ一区二区三区四区 | 国产高清国产精品国产三级| 日韩视频在线欧美| 成人国产av品久久久| 人人妻人人添人人爽欧美一区卜| 久久婷婷青草| 亚洲国产欧美在线一区| videos熟女内射| 亚洲色图 男人天堂 中文字幕 | 插逼视频在线观看| 男男h啪啪无遮挡| 宅男免费午夜| 又粗又硬又长又爽又黄的视频| 天天躁夜夜躁狠狠躁躁| 国产精品国产三级专区第一集| 美国免费a级毛片| 久久久久精品人妻al黑| 亚洲av电影在线观看一区二区三区| 精品亚洲成a人片在线观看| 成人手机av| a级毛片黄视频| 亚洲人成77777在线视频| 久久人人爽av亚洲精品天堂| 午夜视频国产福利| 五月玫瑰六月丁香| av在线老鸭窝| av免费在线看不卡| 精品亚洲乱码少妇综合久久| 成人国产麻豆网| 丰满饥渴人妻一区二区三| 日本黄大片高清| 美女内射精品一级片tv| 女性生殖器流出的白浆| 亚洲,一卡二卡三卡| 赤兔流量卡办理| 亚洲精品色激情综合| 中文字幕亚洲精品专区| 两性夫妻黄色片 | 国产成人精品在线电影| 黄色毛片三级朝国网站| 免费大片黄手机在线观看| 欧美精品av麻豆av| 欧美变态另类bdsm刘玥| 国产一区二区在线观看日韩| 亚洲国产欧美日韩在线播放| 精品福利永久在线观看| 国产精品国产三级国产av玫瑰| 熟女电影av网| 久久久久精品久久久久真实原创| 国产免费视频播放在线视频| 久久久久久久亚洲中文字幕| 日产精品乱码卡一卡2卡三| 卡戴珊不雅视频在线播放| 亚洲欧美色中文字幕在线| 国产精品 国内视频| 人人澡人人妻人| 激情视频va一区二区三区| 宅男免费午夜| 男女啪啪激烈高潮av片| www.熟女人妻精品国产 | 国产伦理片在线播放av一区| 成人手机av| 久久ye,这里只有精品| 一级黄片播放器| 少妇高潮的动态图| 亚洲伊人久久精品综合| 欧美日韩av久久| 亚洲,欧美,日韩| 亚洲av.av天堂| 侵犯人妻中文字幕一二三四区| 五月开心婷婷网| 国产成人精品婷婷| 另类精品久久| 精品熟女少妇av免费看| 中国国产av一级| 日本91视频免费播放| 日韩人妻精品一区2区三区| 伦理电影免费视频| 国产精品久久久久久久久免| 欧美+日韩+精品| 亚洲熟女精品中文字幕| 午夜福利视频精品| 亚洲av在线观看美女高潮| 插逼视频在线观看| 欧美少妇被猛烈插入视频| 中文字幕制服av| 亚洲一码二码三码区别大吗| 日本爱情动作片www.在线观看| 少妇精品久久久久久久| 汤姆久久久久久久影院中文字幕| 美女内射精品一级片tv| 午夜福利乱码中文字幕| 久久精品国产自在天天线| 国产精品久久久久久精品电影小说| 国产成人精品福利久久| 国产欧美另类精品又又久久亚洲欧美| 9热在线视频观看99| 青春草视频在线免费观看| 亚洲国产精品一区二区三区在线| 精品一区在线观看国产| 高清毛片免费看| 日韩av在线免费看完整版不卡| 国产成人a∨麻豆精品| 国产色爽女视频免费观看| av网站免费在线观看视频| √禁漫天堂资源中文www| 日韩 亚洲 欧美在线| 免费看光身美女| 国产精品欧美亚洲77777| 日韩电影二区| 激情五月婷婷亚洲| 国产伦理片在线播放av一区| 日本av手机在线免费观看| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 自拍欧美九色日韩亚洲蝌蚪91| 国产精品麻豆人妻色哟哟久久| 我的女老师完整版在线观看| 亚洲第一av免费看| 大片电影免费在线观看免费| 国产精品麻豆人妻色哟哟久久| 亚洲av男天堂| 久久久a久久爽久久v久久| 最近最新中文字幕免费大全7| 国产亚洲精品第一综合不卡 | 考比视频在线观看| 亚洲av电影在线观看一区二区三区| kizo精华| 久久久久久久精品精品| 女性生殖器流出的白浆| 91精品伊人久久大香线蕉| videosex国产| 国产在视频线精品| 国产成人精品无人区| 大片免费播放器 马上看| 国产精品人妻久久久久久| 欧美日韩精品成人综合77777| 18禁在线无遮挡免费观看视频| 亚洲欧美色中文字幕在线| 22中文网久久字幕| 一级毛片电影观看| 国产无遮挡羞羞视频在线观看| 午夜免费男女啪啪视频观看| 国产一区二区在线观看日韩| 一本大道久久a久久精品| 赤兔流量卡办理| 99香蕉大伊视频| 国产av国产精品国产| 亚洲第一av免费看| 国产精品 国内视频| 少妇猛男粗大的猛烈进出视频| 精品视频人人做人人爽| 亚洲欧美成人综合另类久久久| 久久精品夜色国产| 大陆偷拍与自拍| 这个男人来自地球电影免费观看 | 1024视频免费在线观看| 热re99久久国产66热| 亚洲欧美成人精品一区二区| 日韩大片免费观看网站| 欧美3d第一页| av.在线天堂| 黑人欧美特级aaaaaa片| 男女无遮挡免费网站观看| 2018国产大陆天天弄谢| 午夜视频国产福利| 美女xxoo啪啪120秒动态图| 人人妻人人爽人人添夜夜欢视频| 亚洲 欧美一区二区三区| 九草在线视频观看| videossex国产| 美女脱内裤让男人舔精品视频| av电影中文网址| 狠狠婷婷综合久久久久久88av| 国产精品久久久久成人av| 永久网站在线| 国产黄色视频一区二区在线观看| 久久99精品国语久久久| 亚洲国产欧美日韩在线播放| 国产精品秋霞免费鲁丝片| 伊人久久国产一区二区| 大码成人一级视频| 赤兔流量卡办理| 制服人妻中文乱码| 在线天堂中文资源库| 卡戴珊不雅视频在线播放| 免费高清在线观看日韩| 日本欧美国产在线视频| 日韩av不卡免费在线播放| 国产日韩欧美视频二区| 制服诱惑二区| 老熟女久久久| 一区二区av电影网| 亚洲国产成人一精品久久久| 亚洲综合色惰| 美女中出高潮动态图| 日韩av免费高清视频| 亚洲精华国产精华液的使用体验| 久久人人爽av亚洲精品天堂| 伦理电影大哥的女人| av电影中文网址| 成人毛片a级毛片在线播放| 日韩一区二区三区影片| 国产极品粉嫩免费观看在线| 久久久久久久亚洲中文字幕| av免费在线看不卡| tube8黄色片| 欧美国产精品一级二级三级| 日本黄色日本黄色录像| 国产精品久久久久久av不卡| 中文字幕另类日韩欧美亚洲嫩草| 久久久久国产精品人妻一区二区| 国产黄色免费在线视频| 18禁观看日本| 欧美精品国产亚洲| 深夜精品福利| 成人18禁高潮啪啪吃奶动态图| 久热这里只有精品99| 我的女老师完整版在线观看| 男人爽女人下面视频在线观看| 国产日韩欧美视频二区| 国产探花极品一区二区| 狂野欧美激情性xxxx在线观看| 精品亚洲成a人片在线观看| 两性夫妻黄色片 | 亚洲在久久综合| 9191精品国产免费久久| 母亲3免费完整高清在线观看 | 精品人妻一区二区三区麻豆| 97在线视频观看| 日本黄色日本黄色录像| 成人无遮挡网站| 久久久久久久久久久久大奶| 国产精品蜜桃在线观看| 精品久久蜜臀av无| 午夜福利,免费看| 嫩草影院入口| 国产成人精品无人区| 性色av一级| 中文字幕人妻熟女乱码| 国产又爽黄色视频| 欧美日韩视频高清一区二区三区二| 免费av不卡在线播放| av国产久精品久网站免费入址| 女人久久www免费人成看片| 亚洲中文av在线| 亚洲精品乱久久久久久| 涩涩av久久男人的天堂| 国产一区二区在线观看日韩| 人妻系列 视频| 欧美精品人与动牲交sv欧美| 天堂8中文在线网| 国产色婷婷99| 色婷婷av一区二区三区视频| 黄色视频在线播放观看不卡| 夫妻性生交免费视频一级片| 亚洲欧美清纯卡通| 超色免费av| 亚洲国产精品国产精品| 精品一区在线观看国产| 国产精品一区二区在线观看99| 丝袜脚勾引网站| 婷婷色综合大香蕉| 精品一区二区三卡| av国产精品久久久久影院| 欧美xxⅹ黑人| 国产成人欧美| 五月玫瑰六月丁香| 久久av网站| 精品久久久精品久久久| 蜜臀久久99精品久久宅男| 老熟女久久久| 亚洲欧洲精品一区二区精品久久久 | 国产女主播在线喷水免费视频网站| 黄网站色视频无遮挡免费观看| 精品久久久久久电影网| 国产午夜精品一二区理论片| 卡戴珊不雅视频在线播放| 人成视频在线观看免费观看| 在线亚洲精品国产二区图片欧美| 高清av免费在线| 在线天堂中文资源库| 女性被躁到高潮视频| 波多野结衣一区麻豆| 一级爰片在线观看| 日本av手机在线免费观看| 亚洲欧美精品自产自拍| 香蕉精品网在线| 久久久久视频综合| 久久av网站| 高清欧美精品videossex| 街头女战士在线观看网站| 丝袜喷水一区| 久久久久久久精品精品| 国产女主播在线喷水免费视频网站| 999精品在线视频| 亚洲国产欧美在线一区| 九草在线视频观看| 99久国产av精品国产电影| 亚洲天堂av无毛| 亚洲三级黄色毛片| 久久国产亚洲av麻豆专区| av网站免费在线观看视频| 成年av动漫网址| 亚洲av福利一区| 久久婷婷青草| 日本-黄色视频高清免费观看| 国产精品国产三级国产av玫瑰| 女的被弄到高潮叫床怎么办| 中文字幕人妻熟女乱码| 久久久久久久大尺度免费视频| 男人操女人黄网站| 免费看av在线观看网站| 亚洲精品,欧美精品| 婷婷色综合www| 99视频精品全部免费 在线| 我要看黄色一级片免费的| 欧美日韩一区二区视频在线观看视频在线| av不卡在线播放| 满18在线观看网站| a 毛片基地| 满18在线观看网站| 曰老女人黄片| 亚洲精品乱码久久久久久按摩| 欧美成人午夜免费资源| 久久午夜综合久久蜜桃| 精品熟女少妇av免费看| 亚洲国产成人一精品久久久| 欧美另类一区| 熟女电影av网| 国产精品久久久久久久久免| 色哟哟·www| 成人亚洲欧美一区二区av| 丰满饥渴人妻一区二区三| 亚洲国产欧美在线一区| 伦理电影大哥的女人| 亚洲精品成人av观看孕妇| 亚洲一级一片aⅴ在线观看| 国产在线视频一区二区| 国产精品99久久99久久久不卡 | 亚洲国产精品一区二区三区在线| 18禁观看日本| 在线观看国产h片| 国产精品久久久久久久电影| 晚上一个人看的免费电影| 欧美国产精品一级二级三级| 免费不卡的大黄色大毛片视频在线观看| 久久久久精品人妻al黑| 黄色视频在线播放观看不卡| av片东京热男人的天堂| 免费看光身美女| 狂野欧美激情性bbbbbb| 最近最新中文字幕免费大全7| 国产一区二区三区综合在线观看 | 久久久久精品久久久久真实原创| 免费看不卡的av| 2021少妇久久久久久久久久久| 国产熟女午夜一区二区三区| 在线观看免费视频网站a站| 丰满少妇做爰视频| av免费观看日本| 午夜免费男女啪啪视频观看| 亚洲四区av| av福利片在线| a级片在线免费高清观看视频| 欧美人与性动交α欧美软件 | 久久精品熟女亚洲av麻豆精品| 精品国产国语对白av| 日本色播在线视频| 久久久欧美国产精品| 欧美国产精品va在线观看不卡| 久久这里只有精品19| 亚洲久久久国产精品| 一二三四在线观看免费中文在 |