• <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| 午夜亚洲福利在线播放| 国产区一区二久久| 亚洲av成人一区二区三| 老熟女久久久| 国产av精品麻豆| 国产极品粉嫩免费观看在线| 怎么达到女性高潮| 美女视频免费永久观看网站| 精品乱码久久久久久99久播| 十八禁人妻一区二区| www.精华液| 亚洲一区中文字幕在线| 成年动漫av网址| 少妇裸体淫交视频免费看高清 | 欧美黄色淫秽网站| 亚洲专区字幕在线| 国产淫语在线视频| av超薄肉色丝袜交足视频| 老司机午夜福利在线观看视频| 亚洲午夜理论影院| 精品亚洲成a人片在线观看| 欧美日本中文国产一区发布| 免费在线观看黄色视频的| 最近最新中文字幕大全电影3 | 丝袜美腿诱惑在线| 大码成人一级视频| 午夜激情av网站| 国产又色又爽无遮挡免费看| 色精品久久人妻99蜜桃| 久久影院123| 在线观看免费视频日本深夜| 亚洲久久久国产精品| 桃红色精品国产亚洲av| 国产99久久九九免费精品| 热re99久久精品国产66热6| 一进一出好大好爽视频| 久久九九热精品免费| 成人18禁高潮啪啪吃奶动态图| 一级a爱视频在线免费观看| 亚洲国产看品久久| 久久精品成人免费网站| 色尼玛亚洲综合影院| 一边摸一边抽搐一进一出视频| 国产精品久久久av美女十八| 亚洲专区中文字幕在线| 久久中文字幕人妻熟女| 精品午夜福利视频在线观看一区| 热99国产精品久久久久久7| 村上凉子中文字幕在线| 久久久久久久精品吃奶| 久久影院123| 日韩中文字幕欧美一区二区| 美女高潮到喷水免费观看| 欧美人与性动交α欧美精品济南到| 久久人人爽av亚洲精品天堂| 天天躁夜夜躁狠狠躁躁| 亚洲久久久国产精品| 如日韩欧美国产精品一区二区三区| 两个人看的免费小视频| 飞空精品影院首页| 欧美激情高清一区二区三区| 欧美一级毛片孕妇| 欧美日韩av久久| 精品人妻在线不人妻| av超薄肉色丝袜交足视频| 黑人猛操日本美女一级片| 亚洲欧美精品综合一区二区三区| 中文字幕人妻丝袜一区二区| 欧美成狂野欧美在线观看| 久久人人爽av亚洲精品天堂| 身体一侧抽搐| 国产三级黄色录像| 18禁裸乳无遮挡动漫免费视频| 亚洲自偷自拍图片 自拍| 波多野结衣av一区二区av| 国产亚洲精品久久久久5区| 亚洲成国产人片在线观看| 老司机深夜福利视频在线观看| 国产在线观看jvid| 国产精品久久久久久人妻精品电影| 国产精品免费大片| 妹子高潮喷水视频| 天堂√8在线中文| 久久久国产成人免费| 国产精品电影一区二区三区 | 成人三级做爰电影| 亚洲中文字幕日韩| 午夜福利在线免费观看网站| 大片电影免费在线观看免费| 成人特级黄色片久久久久久久| 亚洲,欧美精品.| 成年人黄色毛片网站| 怎么达到女性高潮| 欧美精品高潮呻吟av久久| 亚洲精品中文字幕一二三四区| 成熟少妇高潮喷水视频| 黄网站色视频无遮挡免费观看| 亚洲精品美女久久久久99蜜臀| 欧美人与性动交α欧美精品济南到| 看片在线看免费视频| av天堂在线播放| 国内毛片毛片毛片毛片毛片| av一本久久久久| 麻豆乱淫一区二区| 91成年电影在线观看| av片东京热男人的天堂| 无限看片的www在线观看| 国产真人三级小视频在线观看| 精品福利永久在线观看| 亚洲色图综合在线观看| 精品久久久久久,| 最近最新中文字幕大全免费视频| 伦理电影免费视频| 国产精品综合久久久久久久免费 | 国产男靠女视频免费网站| а√天堂www在线а√下载 | 精品国产一区二区三区久久久樱花| 欧美丝袜亚洲另类 | 欧美人与性动交α欧美软件| 国产高清激情床上av| 超色免费av| 免费观看精品视频网站| 国产在线观看jvid| 制服诱惑二区| 91九色精品人成在线观看| 国产在线观看jvid| 国产av又大| 在线国产一区二区在线| 亚洲精品久久成人aⅴ小说| 欧美亚洲 丝袜 人妻 在线| 淫妇啪啪啪对白视频| 三上悠亚av全集在线观看| 国产精品av久久久久免费| 成熟少妇高潮喷水视频| 国产高清视频在线播放一区| 男女免费视频国产| 高潮久久久久久久久久久不卡| 大型黄色视频在线免费观看| 亚洲国产欧美一区二区综合| 欧美另类亚洲清纯唯美| 久久精品亚洲精品国产色婷小说| 久久久国产一区二区| 中文字幕av电影在线播放| 亚洲专区字幕在线| 91字幕亚洲| 12—13女人毛片做爰片一| 欧美日韩精品网址| av天堂在线播放| 大陆偷拍与自拍| 亚洲成a人片在线一区二区| 真人做人爱边吃奶动态| 亚洲三区欧美一区| 丰满的人妻完整版| 狠狠婷婷综合久久久久久88av| 亚洲五月婷婷丁香| av网站在线播放免费| 一区二区三区国产精品乱码| 国产亚洲av高清不卡| 午夜影院日韩av| 国产欧美亚洲国产| 久久99一区二区三区| 嫁个100分男人电影在线观看| 成年女人毛片免费观看观看9 | 免费看a级黄色片| 一级片免费观看大全| 精品少妇久久久久久888优播| 韩国精品一区二区三区| 最近最新免费中文字幕在线| 国产高清videossex| 超碰成人久久| 一区二区三区激情视频| 亚洲第一av免费看| 99精品欧美一区二区三区四区| av在线播放免费不卡| 国内毛片毛片毛片毛片毛片| 欧美日韩一级在线毛片| av天堂久久9| 久久青草综合色| 国产99久久九九免费精品| 久久 成人 亚洲| 国产精品国产高清国产av | 夜夜躁狠狠躁天天躁| 久久婷婷成人综合色麻豆| 黄频高清免费视频| 亚洲精品av麻豆狂野| 波多野结衣av一区二区av| 18禁裸乳无遮挡免费网站照片 | 午夜精品久久久久久毛片777| 国产精品亚洲av一区麻豆| 亚洲午夜理论影院| 欧美亚洲日本最大视频资源| 国产真人三级小视频在线观看| 国产高清激情床上av| 国产精品一区二区在线观看99| a级毛片在线看网站| 国产97色在线日韩免费| 黄色a级毛片大全视频| 欧美国产精品va在线观看不卡| 一进一出抽搐gif免费好疼 | 欧美老熟妇乱子伦牲交| 久久中文字幕人妻熟女| 亚洲免费av在线视频| 精品卡一卡二卡四卡免费| 欧美日韩黄片免| 国产欧美日韩一区二区三| 大香蕉久久网| 91国产中文字幕| 久久青草综合色| 国产在线观看jvid| bbb黄色大片| 水蜜桃什么品种好| 午夜福利免费观看在线| 久久人人爽av亚洲精品天堂| 欧美av亚洲av综合av国产av| 精品熟女少妇八av免费久了| 亚洲精品成人av观看孕妇| 国产三级黄色录像| 久久人妻av系列| 婷婷精品国产亚洲av在线 | 身体一侧抽搐| 一边摸一边抽搐一进一小说 | 老鸭窝网址在线观看| 国产1区2区3区精品| 国产主播在线观看一区二区| 黄色片一级片一级黄色片| 激情在线观看视频在线高清 | 亚洲自偷自拍图片 自拍| 99精品欧美一区二区三区四区| 成年女人毛片免费观看观看9 | 18禁裸乳无遮挡免费网站照片 | 午夜福利在线免费观看网站| 久久青草综合色| 精品视频人人做人人爽| 两人在一起打扑克的视频| 国产亚洲精品久久久久久毛片 | 亚洲熟女毛片儿| 俄罗斯特黄特色一大片| 国精品久久久久久国模美| 久久亚洲精品不卡| 狂野欧美激情性xxxx| 老汉色∧v一级毛片| 精品福利永久在线观看| 欧美性长视频在线观看| 亚洲,欧美精品.| 亚洲九九香蕉| 国产xxxxx性猛交| 91成年电影在线观看| x7x7x7水蜜桃| 免费在线观看视频国产中文字幕亚洲| 国产av一区二区精品久久| 久久精品亚洲av国产电影网| 一区二区三区国产精品乱码| av在线播放免费不卡| 欧美乱妇无乱码| 欧美 亚洲 国产 日韩一| 久久狼人影院| 国产精品久久视频播放| 黑人猛操日本美女一级片| 亚洲五月色婷婷综合| 一区二区日韩欧美中文字幕| 亚洲精品中文字幕一二三四区| 亚洲一区高清亚洲精品| 欧美性长视频在线观看| 亚洲欧美色中文字幕在线| 少妇裸体淫交视频免费看高清 | 久久久国产欧美日韩av| 俄罗斯特黄特色一大片| 国产成人免费观看mmmm| 手机成人av网站| 1024香蕉在线观看| 美女 人体艺术 gogo| 99国产精品99久久久久| 午夜福利视频在线观看免费| 日本一区二区免费在线视频| 久久久久久人人人人人| 久久久国产成人免费| 麻豆成人av在线观看| 亚洲七黄色美女视频| 看黄色毛片网站| 十八禁网站免费在线| 国产高清videossex| 欧美在线黄色| 99久久综合精品五月天人人| tube8黄色片| 午夜老司机福利片| 男人舔女人的私密视频| 亚洲一卡2卡3卡4卡5卡精品中文| 久热爱精品视频在线9| 91av网站免费观看| 欧美乱色亚洲激情| 国内毛片毛片毛片毛片毛片| 香蕉久久夜色| 丝袜美腿诱惑在线| 一本综合久久免费| 黄色视频,在线免费观看| 亚洲色图综合在线观看| 99热只有精品国产| 久久热在线av| 深夜精品福利| 国产在线精品亚洲第一网站| 亚洲男人天堂网一区| 欧美日韩亚洲国产一区二区在线观看 | 免费看a级黄色片| 精品熟女少妇八av免费久了| 中文字幕制服av| 制服诱惑二区| 飞空精品影院首页| 午夜福利免费观看在线| 91老司机精品| 亚洲男人天堂网一区| 高清在线国产一区| 国产不卡av网站在线观看| 国产精品免费大片| 两个人免费观看高清视频| 色综合欧美亚洲国产小说| 村上凉子中文字幕在线| 操美女的视频在线观看| 亚洲国产欧美日韩在线播放| 99香蕉大伊视频| 欧美日韩黄片免| 男人的好看免费观看在线视频 | 下体分泌物呈黄色| 夜夜爽天天搞| 99精品在免费线老司机午夜| 国产不卡一卡二| 国产单亲对白刺激| 涩涩av久久男人的天堂| 91国产中文字幕| 亚洲午夜理论影院| 精品久久久精品久久久| 国产精品99久久99久久久不卡| 亚洲av成人av| 女人精品久久久久毛片| 天堂俺去俺来也www色官网| 黄色视频不卡| 中国美女看黄片| 黄色毛片三级朝国网站| 99re在线观看精品视频| 中国美女看黄片| 亚洲av成人一区二区三| 国产精品香港三级国产av潘金莲| 久久精品91无色码中文字幕| 91大片在线观看| 韩国精品一区二区三区| 免费观看人在逋| 丰满人妻熟妇乱又伦精品不卡| 日本五十路高清| 欧美在线一区亚洲| 丰满人妻熟妇乱又伦精品不卡| 午夜精品国产一区二区电影| 91成年电影在线观看| svipshipincom国产片| 妹子高潮喷水视频| 国产有黄有色有爽视频| 成熟少妇高潮喷水视频| 一级黄色大片毛片| 日韩熟女老妇一区二区性免费视频| 成熟少妇高潮喷水视频| 亚洲精品久久成人aⅴ小说| 美女扒开内裤让男人捅视频| 亚洲专区中文字幕在线| 露出奶头的视频| 亚洲,欧美精品.| 成人国语在线视频| 亚洲专区国产一区二区| 国产精品.久久久| 亚洲精品国产一区二区精华液| 18禁观看日本| а√天堂www在线а√下载 | 亚洲国产精品合色在线| 欧美激情极品国产一区二区三区| 大片电影免费在线观看免费| 好男人电影高清在线观看| 国产黄色免费在线视频| 黑人操中国人逼视频| 国产一卡二卡三卡精品| 久久久久久亚洲精品国产蜜桃av| 国产精品亚洲一级av第二区| 大陆偷拍与自拍| 黄色a级毛片大全视频| avwww免费| 视频区图区小说| 国产精品亚洲一级av第二区| 久久人妻福利社区极品人妻图片| 我的亚洲天堂| 80岁老熟妇乱子伦牲交| 日本撒尿小便嘘嘘汇集6| 亚洲精品中文字幕在线视频| 欧美日韩视频精品一区| 水蜜桃什么品种好| 中文字幕人妻丝袜制服| 国内毛片毛片毛片毛片毛片| 欧美人与性动交α欧美精品济南到| 午夜91福利影院| 日韩制服丝袜自拍偷拍| 国产不卡av网站在线观看| 亚洲成a人片在线一区二区| 动漫黄色视频在线观看| 狠狠婷婷综合久久久久久88av| 欧美最黄视频在线播放免费 | 免费在线观看视频国产中文字幕亚洲| 19禁男女啪啪无遮挡网站| 久久久久精品国产欧美久久久| 国产在线一区二区三区精| 青草久久国产| 一二三四在线观看免费中文在| 不卡一级毛片| 岛国毛片在线播放| 精品国产乱子伦一区二区三区| 精品视频人人做人人爽| 久久这里只有精品19| 老汉色∧v一级毛片| 十八禁人妻一区二区| 水蜜桃什么品种好| 大陆偷拍与自拍| 黄色怎么调成土黄色| 亚洲片人在线观看| 亚洲人成电影免费在线| 亚洲精品一卡2卡三卡4卡5卡| 国产有黄有色有爽视频| 日日夜夜操网爽| 欧美激情 高清一区二区三区| 一级片'在线观看视频| 婷婷精品国产亚洲av在线 | 超碰成人久久| 国产精品一区二区在线观看99| 夜夜躁狠狠躁天天躁| 国产欧美日韩一区二区三| 老司机福利观看| 国产日韩欧美亚洲二区| 中文字幕人妻熟女乱码| 一边摸一边做爽爽视频免费| 99国产精品一区二区蜜桃av | videosex国产| 王馨瑶露胸无遮挡在线观看| 欧美在线一区亚洲| 超碰97精品在线观看| 国产成人一区二区三区免费视频网站| 亚洲第一av免费看| 久久人妻av系列| 国产一区二区激情短视频| 91在线观看av| 欧美在线一区亚洲| 欧美日韩黄片免| 女警被强在线播放| 国产一区有黄有色的免费视频| 国产精品欧美亚洲77777| 一本一本久久a久久精品综合妖精| 久久 成人 亚洲| 午夜亚洲福利在线播放| 一二三四社区在线视频社区8| 1024香蕉在线观看| e午夜精品久久久久久久| 日韩制服丝袜自拍偷拍| 香蕉国产在线看| 亚洲国产精品合色在线| 成在线人永久免费视频| 精品少妇久久久久久888优播| 国产成人系列免费观看| 中文字幕人妻熟女乱码| 久久久精品区二区三区| 在线观看午夜福利视频| 少妇猛男粗大的猛烈进出视频| 成年人午夜在线观看视频| 成在线人永久免费视频| 人妻丰满熟妇av一区二区三区 | 两个人免费观看高清视频| 亚洲欧美激情综合另类| 国产亚洲一区二区精品| 精品福利永久在线观看| 热99国产精品久久久久久7| 久久国产精品大桥未久av| 人妻丰满熟妇av一区二区三区 | 免费观看a级毛片全部| 日本黄色视频三级网站网址 | 在线观看免费视频日本深夜| 亚洲国产欧美一区二区综合| 国产精品国产高清国产av | 高清欧美精品videossex| 黄色视频,在线免费观看| 午夜91福利影院| 日韩视频一区二区在线观看| 在线国产一区二区在线| 久久久国产精品麻豆| 一级毛片高清免费大全| 在线天堂中文资源库| 黄色怎么调成土黄色| 欧美黑人欧美精品刺激| 咕卡用的链子| 免费人成视频x8x8入口观看| 欧美日韩亚洲高清精品| 久久久久国产一级毛片高清牌| 午夜免费观看网址| 久久精品国产综合久久久| 亚洲五月天丁香| 天天躁夜夜躁狠狠躁躁| www.精华液| 黄色视频不卡| 老熟妇乱子伦视频在线观看| 欧美日韩一级在线毛片| 国产精品香港三级国产av潘金莲| 国产99久久九九免费精品| 日本a在线网址| 激情视频va一区二区三区| x7x7x7水蜜桃| 亚洲视频免费观看视频| 欧美日韩国产mv在线观看视频| 一级,二级,三级黄色视频| 少妇粗大呻吟视频| 欧美另类亚洲清纯唯美| 大香蕉久久网| 欧美精品亚洲一区二区| 久久久精品国产亚洲av高清涩受| 久久久精品免费免费高清| 91精品国产国语对白视频| 精品福利永久在线观看| 大型黄色视频在线免费观看| 国产精品 欧美亚洲| 欧美日韩中文字幕国产精品一区二区三区 | 天堂中文最新版在线下载| 水蜜桃什么品种好| 一进一出抽搐动态| 日韩精品免费视频一区二区三区| 啦啦啦视频在线资源免费观看| 亚洲av日韩在线播放| 悠悠久久av| 久久久国产一区二区| 国产午夜精品久久久久久| 国产精品乱码一区二三区的特点 | 一级,二级,三级黄色视频| 好男人电影高清在线观看| svipshipincom国产片| 久久久国产精品麻豆| 亚洲专区国产一区二区| 日韩一卡2卡3卡4卡2021年| 亚洲视频免费观看视频| 香蕉久久夜色| 午夜精品在线福利| 亚洲国产看品久久| 一级毛片高清免费大全| 久久精品成人免费网站| 午夜影院日韩av| 欧美久久黑人一区二区| 国产伦人伦偷精品视频| 久久久久久免费高清国产稀缺| 中文字幕色久视频| 日韩三级视频一区二区三区| 成年人午夜在线观看视频| 操出白浆在线播放| av视频免费观看在线观看| 黑人操中国人逼视频| 欧美日韩成人在线一区二区| 国产一区有黄有色的免费视频| 人妻丰满熟妇av一区二区三区 | 精品久久久精品久久久| 欧美乱妇无乱码| 美女午夜性视频免费| 国产成人免费无遮挡视频| 麻豆av在线久日| 丝袜在线中文字幕| 久久香蕉激情| 91字幕亚洲| 欧美日韩亚洲综合一区二区三区_| 大码成人一级视频| 日韩 欧美 亚洲 中文字幕| 777米奇影视久久| 欧美黑人精品巨大| xxx96com| 五月开心婷婷网| 欧美人与性动交α欧美软件| 久久亚洲真实| 国产无遮挡羞羞视频在线观看| 日韩有码中文字幕| tube8黄色片| 高清在线国产一区| 精品人妻熟女毛片av久久网站| 中文字幕高清在线视频| 日本a在线网址| 国产视频一区二区在线看| 99精品久久久久人妻精品| 精品亚洲成国产av| 欧洲精品卡2卡3卡4卡5卡区| 午夜免费观看网址| 丁香六月欧美| 狠狠婷婷综合久久久久久88av| 国产精品98久久久久久宅男小说| 香蕉国产在线看| 精品人妻1区二区| 人妻久久中文字幕网| 久久天堂一区二区三区四区| 99精品在免费线老司机午夜| 精品国产乱码久久久久久男人| 欧美精品亚洲一区二区| 后天国语完整版免费观看| 99久久人妻综合| 欧美中文综合在线视频| 国产一区二区三区在线臀色熟女 | 亚洲国产毛片av蜜桃av| 啦啦啦在线免费观看视频4| 视频在线观看一区二区三区| 亚洲国产精品sss在线观看 | 下体分泌物呈黄色| 操美女的视频在线观看| 久久久水蜜桃国产精品网| 黑人巨大精品欧美一区二区mp4| 十分钟在线观看高清视频www| 19禁男女啪啪无遮挡网站| 久久国产精品影院| 91麻豆av在线| 天天添夜夜摸| 老司机靠b影院|