江葉楓, 孫 凱, 郭 熙*, 葉英聰, 饒 磊, 李偉峰

1.江西農(nóng)業(yè)大學國土資源與環(huán)境學院, 江西 南昌 330045 2.江西省鄱陽湖流域農(nóng)業(yè)資源與生態(tài)重點實驗室, 江西 南昌 330045

基于環(huán)境因子和鄰近信息的土壤屬性空間分布預測

江葉楓1,2, 孫 凱1,2, 郭 熙1,2*, 葉英聰2, 饒 磊1,2, 李偉峰1,2

1.江西農(nóng)業(yè)大學國土資源與環(huán)境學院, 江西 南昌 330045 2.江西省鄱陽湖流域農(nóng)業(yè)資源與生態(tài)重點實驗室, 江西 南昌 330045

為探索鄉(xiāng)鎮(zhèn)尺度上土壤屬性空間分布預測的最佳方法，以江西省萬年縣齊埠鎮(zhèn)為例，借助四方位搜索法、地統(tǒng)計學和遙感影像分析技術(shù)提取環(huán)境因子(地形因子和植被覆蓋指數(shù))和鄰近信息〔w(有機質(zhì))與w(速效鉀)〕，構(gòu)建OK法(普通克里金法)、RK1法(僅基于環(huán)境因子的回歸克里金法)以及RK2法(基于環(huán)境因子和鄰近信息的回歸克里金法)對齊埠鎮(zhèn)耕地表層(0～20 cm)土壤w(有機質(zhì))、w(速效鉀)空間分布進行預測. 結(jié)果表明：齊埠鎮(zhèn)土壤w(有機質(zhì))平均值為35.03 g/kg，w(速效鉀)平均值為96.73 mg/kg，均為中等空間變異性. 對62個樣點進行建模，16個測試樣點進行獨立驗證的誤差分析表明，RK2法對土壤w(有機質(zhì))、w(速效鉀)預測結(jié)果的均方根誤差、平均絕對誤差和平均相對誤差較OK法分別降低了18.05%、18.01%、21.77%和7.25%、9.49%、9.84%；較RK1法分別降低了22.48%、20.91%、22.02%和9.27%、12.61%、13.52%. 研究顯示，RK2法明顯提高了土壤w(有機質(zhì))、w(速效鉀)空間分布模擬精度，并且存在改進和提高的空間.

土壤屬性； 四方位搜索法； 回歸克里金法； 空間預測

Abstract： As the most important determinants of soil quality, soil properties significantly influence farmland use and soil environmental protection. In order to find the best interpolation method of soil properties, including soil organic matter (SOM) and Soil Available Potassium (AK) in Qibu Town in Wannian County, Jiangxi Province, Ordinary Kriging (OK), Regression Kriging 1 (RK1, based on environmental factors), Regression Kriging 2 (RK2, based on environmental factors and neighbor information) were adopted to predict the distribution of soil properties. Environmental factors were extracted by digital terrain and remote sensing image analysis technique. The four-direction search method was applied to get the neighbor information. To establish and validate the three methods, 78 points of surface soil samples (0-20 cm) were collected and randomly divided into two groups, as modeling points (62) and validation points (16). The results showed that SOM content ranged from 17.30-53.58 g/kg, with an average value of 35.03 g/kg, a moderate variability. The nugget/sill ratio was 0.59, indicating a moderate spatial dependence for SOM. AK content ranged from 50.98 to 152.13 mg/kg, with an average value of 96.73 mg/kg, a moderate variability. The nugget/sill ratio was 0.66, indicating a moderate spatial dependence for AK. The prediction map obtained by RK2 model was more consistent with the true geographical information than OK and RK1. Moreover, RK2 model reduced the prediction errors. Compared to RK1 model, the root mean square errors, the mean absolute errors and the mean relative errors of RK2 were 18.05%, 18.01%, and 21.77% (for SOM), and 7.25%, 9.49%, and 9.84% (for AK), smaller than those of OK. The same numbers for the OK model with the validation points were 22.48%, 20.91%, and 22.02% (for SOM), and 9.27%, 12.61%, and 13.52% (for AK), smaller than those of RK1. The results suggested that it is helpful for improving the prediction accuracy to employ neighbor information and environmental factors in spatial prediction of soil properties. Therefore, RK2 could be recognized as the best interpolation method, but could be improved in the future.

Keywords： soil properties; four-direction search method; Regression Kriging method; spatial prediction

土壤屬性是影響土壤質(zhì)量的最重要因素，對農(nóng)田利用和土壤環(huán)境保護具有重要影響. 受生物[1]、氣候[2]、母質(zhì)[3]和地形[4]等因素影響導致其在空間分布上呈現(xiàn)非均勻性. 快速準確的獲取土壤屬性空間分布是有效管理土壤、科學規(guī)劃土地利用、精準農(nóng)業(yè)施肥的基礎(chǔ)[5]. 目前基于實地采樣獲取的土壤屬性信息遠遠不能滿足精準農(nóng)業(yè)和環(huán)境保護的實際需求. 由于成本限制大規(guī)模采樣并不現(xiàn)實，且仍然無法得到研究區(qū)域空間上土壤屬性含量. 因此，預測土壤屬性空間分布一直是土壤學研究的熱點問題之一[6-8].

目前，針對土壤屬性預測方法有很多，其中OK法(Ordinary Kriging，OK)基于鄰近相關(guān)采樣點權(quán)重來對未知點進行預測[9]，而且能夠?qū)︻A測誤差進行理論估計，是應用最普遍也最廣泛的方法[10]. 但土壤屬性的涵養(yǎng)和運移分布受到地形、植被等環(huán)境因子的影響，具有高度的空間變異性[11]. Eldeiry等[12]研究表明，僅基于鄰近相關(guān)采樣點權(quán)重而不考慮地形、植被覆蓋指數(shù)等環(huán)境因子模擬土壤屬性空間分布的OK法預測精度已經(jīng)受到了較大限制. 為此，學者們引入RK法(Regression Kriging, RK)將環(huán)境因子作為輔助變量納入到土壤屬性空間分布預測中. 如張素梅等[13-14]利用RK法將地形和植被覆蓋指數(shù)等環(huán)境因子納入模型中預測土壤屬性空間分布并取得了較好的結(jié)果. 但是，由于以往研究成果都集中在縣域尺度上[15-18]，而鄉(xiāng)鎮(zhèn)尺度上環(huán)境因子相對單一，與縣域尺度相比較模糊；且土壤屬性具有高度的空間變異性，僅基于環(huán)境因子的RK法這一類全局模型無法捕捉到土壤屬性變異的局部特征，擬合系數(shù)較低，預測精度再次受到限制. 近年來，越來越多的學者將鄰近信息納入土壤屬性的空間預測中來. 徐劍波等[19]利用歐式距離來尋找離建模點最近的5個土壤屬性值，在不同樣本下的誤差較OK法都有明顯的降低. 李啟權(quán)等[20]根據(jù)空間距離衰減規(guī)律，選取離插值點最近的3個土壤樣點的屬性值，在不同土壤屬性間較OK法誤差降低達到顯著水平. 僅基于環(huán)境因子進行逐步多元回歸只能在全局范圍內(nèi)進行預測[21]；而通過把鄰近信息納入土壤屬性的空間分布預測，能夠更好地揭示土壤屬性空間分布的局部變化特征，能在提高精度的同時更加真實的反映土壤屬性空間變異情況.

為此，該研究以江西省萬年縣齊埠鎮(zhèn)耕地土壤w(有機質(zhì))、w(速效鉀)為研究對象. 根據(jù)反距離加權(quán)算法思想[22]和空間自相關(guān)理論[23]，運用四方位搜索法尋找離插值點最近的w(有機質(zhì))、w(速效鉀)作為鄰近信息(P1、P2、P3、P4)；同時利用DEM(數(shù)字高程模型)數(shù)據(jù)與遙感影像分析技術(shù)提取地形、植被覆蓋指數(shù)等環(huán)境因子；構(gòu)建RK2法預測齊埠鎮(zhèn)土壤w(有機質(zhì))、w(速效鉀)空間分布，并與RK1(僅基于環(huán)境因子的回歸克里金法)法和OK法進行比較. 以期為鄉(xiāng)鎮(zhèn)尺度上土壤屬性空間分布預測提供方法參考.

1 材料與方法

1.1 研究區(qū)域概況

該研究以江西省萬年縣齊埠鎮(zhèn)為研究區(qū). 該區(qū)位于萬年河下游，地處余干、萬年和鄱陽三縣交匯之處，是萬年縣的西大門，屬于亞熱帶季風性氣候，總面積為72.8 km2. 境內(nèi)水資源比較豐富，年均降水量為 1 766 mm，無霜期為263 d. 地貌類型以平原為主. 土地利用類型主要為耕地和林地，其中耕地面積23.4 km2，占總面積的32.17%；林地面積40.8 km2，占總面積的56.04%. 成土母質(zhì)有第四紀紅色黏土、碳酸鹽類風化物、泥頁巖類風化物. 土壤類型主要是紅壤.

1.2 土壤采樣

土壤樣品采集于2014年8月. 按照均勻性、代表性和連續(xù)性的原則，同時考慮地形部位的基礎(chǔ)上，采用“S”形采樣、多點混合的方法采集0～20 cm的耕地表層土壤樣品，經(jīng)過充分混合最后留取1 kg土樣，得到78個耕地表層土壤樣品(其中水田69個，旱地9個). 土壤樣品經(jīng)過自然風干后，帶回實驗室磨碎過篩，采用重鉻酸鉀(K2Cr2O7)油浴加熱法測定土壤w(有機質(zhì))；采用乙酸銨浸提-火焰光度法測定土壤w(速效鉀)[24]. 為驗證筆者所用方法對研究區(qū)土壤屬性的空間分布預測精度，隨機均勻選取62個采樣點作為建模點，用于插值；剩下16個作為測試點，不參與預測過程只在后續(xù)精度分析中用于驗證預測精度.

1.3 環(huán)境因子提取

土壤屬性的空間分布是物理、化學和生物的綜合作用結(jié)果[13]. 在中小尺度上，主要受定性和定量變量的影響[15,25]，其中定量變量主要有地形因子和植被覆蓋指數(shù). 參考國內(nèi)外的研究成果[26-28]，該研究選取高程(H)、坡度(S)、坡向(AS)、曲率(C)、坡度變率(SOS)、坡向變率(SOA)、地形起伏度(QFD)、河流動能指數(shù)(Ω)、地形濕度指數(shù)(TI)和植被覆蓋指數(shù)(NDVI)10個因素作為影響土壤屬性空間分布的環(huán)境因子(見圖1). 其中植被覆蓋指數(shù)由GF-1(拍攝日期為2015年11月15日，空間分辨率8 m)的第三波段和第四波段在ArcGIS 10.2中進行柵格計算獲??；其他地形因子可通過GIS空間分析從DEM數(shù)據(jù)(由1∶10 000 地形圖生成)獲得. 各環(huán)境因子的計算公式見文獻[13].

圖1 研究區(qū)內(nèi)影響土壤屬性空間分布的主要環(huán)境因子及采樣點分布Fig.1 The main environmental factors affecting the spatial distribution of soil properties and the distribution of sampling sites in the study area

1.4 鄰近信息提取

利用四方位搜索法(the four-direction search method)獲取插值點附近土壤屬性值〔w(有機質(zhì))、w(速效鉀)〕. 四方位搜索法是基于反距離加權(quán)算法思想和空間自相關(guān)理論，根據(jù)采樣點之間的距離越近影響越大，距離越遠影響越小的思想，以采樣點為坐標原點，變程a為半徑，在每個象限內(nèi)選擇一個鄰近點土壤屬性值P1、P2、P3、P4. 四方位搜索法步驟：

a) 在MATLAB中新建一個mat文件，導入所有采樣點地理坐標(x，y)和土壤屬性值(z)；

b) 從文件中隨機抽取一個點，賦值給Mi；

c) 在mat文件剩余點中隨機抽取一個點，賦值給Ni；

d) 計算M、N兩點之間的歐氏距離，設為d，令Temp1=Temp2=Temp3=Temp4=a，若d>a，則回到步驟3)；

e) 若Ni(x)>Mi(x),Ni(y)≥Mi(y)且d≤Temp1，則Temp1=d,P1=Ni(z)；

f) 若Ni(x)≤Mi(x),Ni(y)>Mi(y)且d≤Temp2，則Temp2=d,P2=Ni(z)；

g) 若Ni(x)

h) 若Ni(x)≥Mi(x),Ni(y)

i) 重復步驟b)～h)，直至遍歷mat文件中所有點，將P1、P2、P3、P4存放在文本中.

四方位搜索法在MATLABR2014a中通過編程實現(xiàn)，出現(xiàn)的空值由其他象限的平均值進行計算. 其中P1、P2、P3、P4分別代表以采樣點為坐標原點，4個象限內(nèi)離采樣點最近點的土壤w(有機質(zhì))或w(速效鉀)；Mi和Ni表示兩個采樣點的地理坐標，Mi(x)和Ni(x)表示兩個采樣點的經(jīng)度，Mi(y)和Ni(y)表示兩個采樣點的緯度.

1.5 研究方法

1.5.1 RK2法

在鄉(xiāng)鎮(zhèn)尺度上，RK1法往往精度不理想. 研究[19-20]表明，加入土壤屬性的鄰近信息后可以提高空間分布模擬精度. 因此，要較為精確地模擬土壤屬性的空間分布情況，必須綜合考慮土壤屬性的空間自相關(guān)性和影響土壤屬性空間變異的環(huán)境因子. RK2法是通過多元逐步回歸建立起土壤屬性和環(huán)境因子、鄰近信息之間的回歸預測方程，根據(jù)回歸方程擬合出一個代表確定性部分的趨勢項表面，同時對回歸預測殘差進行OK法插值得到代表隨機性部分的殘差項，最后將二者進行疊加. 其過程可表示為

式中，f1(x)為RK2法在x處的模擬值，g1(x)為RK2法在x處根據(jù)多元逐步回歸方程擬合的趨勢項，m1(x) 為OK法在x處插值的殘差項.

RK2法預測過程：首先在SAS軟件中應用多元逐步回歸得到土壤w(有機質(zhì))、w(速效鉀)與環(huán)境因子以及鄰近信息的回歸方程，利用該方程得到土壤w(有機質(zhì))、w(速效鉀)的趨勢項表面；運用GS+(version 7)對回殘差進行半方差分析，得到殘差的半方差模型，然后運用OK法插值得到代表隨機性部分的殘差項. 最后將趨勢項和殘差項在ArcGIS 10.2中進行柵格疊加，得到基于環(huán)境因子和鄰近信息的土壤屬性空間分布預測結(jié)果.

1.5.2 參照方法

將OK法和RK1法作為參照方法. OK法作為傳統(tǒng)的地統(tǒng)計學方法被廣泛的應用于土壤屬性的空間插值中，是單個變量的局部線性最優(yōu)無偏估計方法；也是應用最廣泛、最穩(wěn)健也最常用的一種方法. RK1法與RK2法步驟一樣，僅僅是在進行逐步回歸時只將環(huán)境因子與土壤屬性進行回歸分析，而不把鄰近信息納入逐步回歸方程.

1.6 預測精度分析

通過ArcGIS 10.2得到建模集和測試集，以RMSE(均方根誤差)、MAE(平均絕對誤差)、MRE(平均相對誤差)對建模和測試的土壤w(有機質(zhì))、w(速效鉀) 預測值與實際采樣值進行對比分析，得出精度評價結(jié)果，其公式分別為

(1)

(2)

(3)

2 結(jié)果與討論

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

研究區(qū)土壤屬性的描述性統(tǒng)計特征分析見表1. 研究區(qū)土壤w(有機質(zhì))的平均值為35.03g/kg，值域范圍17.30～53.58g/kg；w(速效鉀)的平均值為96.73mg/kg，值域范圍50.98～152.13mg/kg. 從K-S檢驗可以判斷兩種土壤屬性統(tǒng)計值均符合正態(tài)分布. 從變異系數(shù)〔w(有機質(zhì))，23.61%；w(速效鉀)，21.72%〕可以看出，均為中等變異性.

表1 研究區(qū)土壤屬性的描述性統(tǒng)計特征

數(shù)據(jù)的正態(tài)性是OK法的前提，只有當數(shù)據(jù)服從正態(tài)分布時，OK法才有效[10]. 由圖2可見，RK1法和RK2法預測土壤w(有機質(zhì))、w(速效鉀)的殘差值基本沿著正態(tài)線分布，符合正態(tài)分布，可以進行半方差分析和OK法插值.

2.2 土壤屬性與環(huán)境因子及鄰近信息的相關(guān)性

研究區(qū)土壤屬性與環(huán)境因子及鄰近信息的相關(guān)性如表2所示. 由表2可見，土壤w(有機質(zhì))與地形濕度指數(shù)呈顯著的負相關(guān)，相關(guān)系數(shù)為-0.229(P<0.05)，表明地形低洼處土壤w(有機質(zhì))越低；與其他環(huán)境因子沒有顯著相關(guān)性. 這與以往研究結(jié)果[13,16]不一致，究其原因：①現(xiàn)有研究成果都集中在縣域尺度且地形多為丘陵和山地，地貌類型多樣，溝壑縱橫，環(huán)境差異大；而在鄉(xiāng)鎮(zhèn)尺度和平原上可能因為范圍小、地勢平坦導致由環(huán)境等宏觀因素引起的土壤w(有機質(zhì))的涵養(yǎng)與運移過程被削弱[29]；②該研究采樣點較少，無法比較在更多樣點分布下的分析結(jié)果，這也在一定程度上減弱了環(huán)境因子與土壤w(有機質(zhì))的相關(guān)性[30].w(速效鉀)與坡度、坡度變率呈顯著負相關(guān)，相關(guān)系數(shù)分別為-0.256(P<0.05)和-0.266(P<0.05)，表明坡度越陡，坡度變化越大，w(速效鉀)越低，這與以往研究結(jié)果[5]一致. 地形濕度指數(shù)和曲率與w(速效鉀)的相關(guān)性均達到顯著(R=0.301，P<0.05；R=0.263，P<0.05). 其他地形因子與w(速效鉀)的相關(guān)性未達到顯著水平. 此外，鄰近信息P1、P2、P3、P4與土壤屬性均有極顯著的相關(guān)性(P<0.01)，相關(guān)系數(shù)均達到0.4以上，表明在鄉(xiāng)鎮(zhèn)尺度上土壤屬性存在較強的自相關(guān)性，鄰近信息對土壤屬性有著非常重要的影響.

圖2 土壤w(有機質(zhì))、w(速效鉀)RK1法和RK2法預測殘差的Q-Q圖Fig.2 Q-Q plot of prediction residual by RK1 and RK2 for soil organic matter and Available Potassium

影響因素w(有機質(zhì))w(速效鉀)影響因素w(有機質(zhì))w(速效鉀)植被覆蓋指數(shù)0.0170.078地形起伏度-0.045-0.208高程-0.149-0.192河流動能指數(shù)-0.028-0.092坡度-0.085-0.256地形濕度指數(shù)-0.229*0.301*坡向-0.045-0.137P10.674**0.623**曲率0.0620.263*P20.430**0.533**坡度變率0.065-0.266*P30.482**0.523**坡向變率0.1150.059P40.551**0.421**

注：*表示P<0.05；**表示P<0.01.

2.3 半方差函數(shù)分析

在GS+(version 7)中對土壤屬性、RK1法殘差和RK2法殘差進行半方差函數(shù)的擬合，用半方差函數(shù)描述其空間變異性. 由表3可見，除w(有機質(zhì))和w(速效鉀)RK1法預測殘差的理論模型為指數(shù)外，其他數(shù)據(jù)項最優(yōu)模型均為球狀模型，各模型的擬合系數(shù)均達到0.8以上. 從模型的參數(shù)來看，土壤w(有機質(zhì))、w(速效鉀)的塊金效應值均小于0.75，變程分別為3 503和1 791 m，表明二者均具有較強的自相關(guān)性，空間自相關(guān)的范圍較大. 4種回歸預測殘差的半方差模型參數(shù)與原變量變化較小，基本保留了原變量的空間結(jié)構(gòu)特征.

2.4 土壤屬性的多元逐步回歸模擬過程

多元逐步回歸能保證每個鄰近信息和環(huán)境因子進入回歸模型且能去除自變量間的共線性[21,28]. 考慮到土壤屬性空間分布的位置關(guān)系，在逐步回歸擬合時把地理坐標(x，y)也作為自變量進入回歸方程. 表4是土壤屬性的多元逐步回歸模擬過程. 在土壤w(有機質(zhì))的RK1法預測中，地形濕度指數(shù)、坡向變率兩個環(huán)境因子和地理坐標x參與了土壤w(有機質(zhì))的空間分布模擬；在RK2法預測中坡向變率和鄰近信息P1參與了土壤w(有機質(zhì))的空間分布模擬. 在土壤w(速效鉀)RK1空間分布中，植被覆蓋指數(shù)、地形濕度指數(shù)、河流動能指數(shù)3個環(huán)境因子參與了速效鉀的空間分布模擬；在RK2中僅鄰近信息P1參與了速效鉀的空間分布模擬.

表3 土壤屬性、RK1法殘差和RK2法殘差的半方差函數(shù)參數(shù)

表4 研究區(qū)土壤屬性含量的多元逐步回歸擬合過程

從各模型的擬合系數(shù)來看，土壤w(有機質(zhì))、w(速效鉀)RK1法的擬合系數(shù)偏低，這與已有研究結(jié)果[5,8]相似. 將OK法、RK1法和RK2法對土壤w(有機質(zhì))、w(速效鉀)的預測結(jié)果與實測值進行相關(guān)性分析，結(jié)果表明RK2法對土壤w(有機質(zhì))、w(速效鉀)的預測值與真實值的相關(guān)系數(shù)要大于OK法和RK1法(見表5). 其中，RK2法對土壤w(有機質(zhì))的建模點和測試點的預測結(jié)果與采樣點實測值的相關(guān)系數(shù)較RK1法分別提高了0.164和0.244，較OK法分別提高了0.013和0.259. RK2法對w(速效鉀)的建模點和測試點的預測結(jié)果與采樣點實測值的相關(guān)系數(shù)較RK1法分別提高了0.089和0.169，較OK法分別提高了0.052和0.201. 這表明在鄉(xiāng)鎮(zhèn)尺度下，加入鄰近信息的RK2法較僅基于環(huán)境因子的RK1法和OK法能夠更加真實的反映土壤屬性的空間分布，能夠有效的揭示土壤屬性的復雜空間變異關(guān)系(見圖3). 而RK1法對w(有機質(zhì))建模點的預測結(jié)果與采樣點實測值的相關(guān)系數(shù)較OK法降低了0.151；RK1法對w(速效鉀)建模點的預測結(jié)果與采樣點實測值的相關(guān)系數(shù)較OK法降低了0.037. 這是因為回歸克里金的預測精度取決于回歸預測方程的決定系數(shù)，而在鄉(xiāng)鎮(zhèn)等小尺度上，環(huán)境因素不如縣域那么明顯[25]，導致僅由環(huán)境因子進行回歸建模的RK1法預測結(jié)果與實測值得相關(guān)系數(shù)偏低(見表4)，甚至低于基于土壤屬性空間自相關(guān)程度(塊金效應)[31-32]的OK法. 因此，在鄉(xiāng)鎮(zhèn)尺度上有必要將鄰近信息納入土壤屬性空間分布預測中.

2.5 預測精度評價

表5 不同預測方法預測結(jié)果與實測值的相關(guān)性

注：*表示P<0.05；**表示P<0.01.

土壤屬性預測精度分析結(jié)果(見表6)表明，引入鄰近信息進行土壤屬性空間分布預測的RK2法明顯優(yōu)于僅基于環(huán)境因子進行空間分布預測的RK1法和僅基于采樣點進行空間內(nèi)插的OK法. 從建模集預測效果看，與OK法進行比較，RK2法對土壤w(有機質(zhì))預測結(jié)果的RMSE、MAE、MRE分別降低了7.28%、9.88%、13.39%；與RK1法進行比較，分別降低了21.15%、23.06%、32.07%. 與OK法進行比較，RK2法對w(速效鉀)預測結(jié)果的RMSE、MAE、MRE分別降低了5.90%、6.80%、7.11%；與RK1法進行比較，分別降低了8.17%、8.21%、9.99%. 從測試集預測效果看，與OK法進行比較，RK2法對土壤w(有機質(zhì))預測結(jié)果的RMSE、MAE、MRE分別降低了18.05%、18.01%、21.77%；與RK1法進行比較，分別降低了22.48%、20.91%、22.02%. 與OK法進行比較，RK2法對w(速效鉀)預測結(jié)果的RMSE、MAE、MRE分別降低了7.25%、9.49%、9.84%；與RK1法進行比較，分別降低了9.27%、12.61%、13.52%. 基于環(huán)境因子和鄰近信息的方法預測精度得到了較為明顯的提高.

表6 土壤屬性預測方法精度對比

圖3 不同回歸方法預測的土壤w(有機質(zhì))和w(速效鉀)的空間分布模擬結(jié)果Fig.3 Maps of soil organic matter and available potassium by different regression methods

將OK法、RK1法和RK2法的預測空間分布結(jié)果以相同分辨率在ArcGIS 10.2中顯示(見圖3)，3種方法在土壤屬性空間分布預測中差異顯著. OK法預測土壤屬性的空間分布呈現(xiàn)“南北低、中部高”的總體趨勢，這與DEM數(shù)據(jù)變化情況較吻合. 從空間分布模擬圖來看，土壤屬性高低值分布界限較RK1法相對分明，而RK1法得到的空間分布模擬結(jié)果較平滑，高低值界限不清晰，難以體現(xiàn)土壤屬性的空間變異性. RK2法得到空間模擬圖的高低值呈塊狀分布，體現(xiàn)了土壤屬性的空間變異信息；同時研究區(qū)內(nèi)出現(xiàn)了較多高值區(qū)域包含的低值部分，結(jié)合圖1和圖3可知該低值區(qū)域DEM要高于周圍區(qū)域，這可能是因為齊埠鎮(zhèn)汛期多降雨導致高程較高處受暴雨沖刷相對嚴重，從而使得土壤w(有機質(zhì))、w(速效鉀)流失，這比較符合土壤屬性地學運動規(guī)律與實際情況，揭示了土壤屬性隨地形變化的細節(jié)信息. 從預測結(jié)果的值域范圍來看，RK2法預測值域范圍最為接近統(tǒng)計分析值〔w(有機質(zhì))為21.83～48.53 g/kg，w(速效鉀)為63.75～123.65 mg/kg，見圖3〕，RK1法預測值域范圍最小〔w(有機質(zhì))為26.50～42.62 g/kg，w(速效鉀) 為64.23～116.30 mg/kg，見圖3〕. 該研究表明僅基于環(huán)境因子的預測結(jié)果更趨向于平均值，難以克服平滑效應；而同時加入鄰近信息和環(huán)境因子的RK2法能夠更為精確的預測土壤屬性的空間分布. 土壤屬性空間變異為突變而非漸變[33]，在鄉(xiāng)鎮(zhèn)尺度上，土壤屬性的空間自相關(guān)性范圍大(見表3)，同時環(huán)境因子也不如縣域那么明顯[25]，RK1法和OK法難以表達土壤屬性含量突變的區(qū)域. 而RK2法根據(jù)土壤屬性的空間自相關(guān)性和影響土壤屬性的宏觀因素把鄰近信息和環(huán)境因子同時加入回歸方程，能夠有效的提高土壤屬性空間分布的預測精度，使土壤屬性空間分布模擬更加符合實際情況.

2.6 討論

土壤屬性的空間分布受眾多因素的綜合作用. 該研究結(jié)果顯示，在鄉(xiāng)鎮(zhèn)尺度上受環(huán)境因子和鄰近信息的綜合作用，但鄰近信息影響更大. 相關(guān)性分析結(jié)果表示鄰近信息與土壤屬性含量的相關(guān)系數(shù)均大于環(huán)境因子. 這是因為在鄉(xiāng)鎮(zhèn)等小尺度下土壤屬性空間自相關(guān)性較強，空間連續(xù)性范圍大；這與以往研究結(jié)果一致. 而大部分環(huán)境因子在鄉(xiāng)鎮(zhèn)尺度上與土壤屬性含量相關(guān)性并不顯著，這與張素梅等[13,15]研究結(jié)果不一致，究其原因：①齊埠鎮(zhèn)高程均<200 m，屬于典型的平原地帶，地勢起伏相對平緩；同時由于齊埠鎮(zhèn)森林覆蓋率達到56%以上，地表團聚性強，土壤w(有機質(zhì))、w(速效鉀)等屬性的涵養(yǎng)和運移過程及分布受環(huán)境因子變化影響較小. ②齊埠鎮(zhèn)耕地圖斑分布相對分散，目前采樣的方式和數(shù)量難以反映土壤屬性的真實分布情況，其他隨機性因素如耕作、施肥、秸稈還田等在目前的采樣密度內(nèi)對土壤屬性含量影響較大，后續(xù)應加大采樣密度和改進采樣方式才能更加真實的反映鄉(xiāng)鎮(zhèn)尺度上土壤屬性含量與環(huán)境因子的相關(guān)性.

在現(xiàn)有土壤屬性空間分布預測中，OK法是應用最普遍也最具代表性的，雖然能夠通過鄰近相關(guān)采樣點土壤屬性值權(quán)重來預測未知點、對預測結(jié)果給出誤差等優(yōu)點，但該方法沒有考慮環(huán)境因子的影響，導致其制圖效果和預測精度上還不是很理想. RK1法和RK2法不僅能夠構(gòu)建土壤屬性與自變量之間的回歸預測方程，通過回歸預測方程擬合出代表確定性部分的趨勢項，而且在殘差存在空間自相關(guān)的情況下充分利用RK1和RK2殘差進行OK法插值，考慮了采樣點的結(jié)構(gòu)性和隨機性的空間模擬. RK1法僅基于影響土壤屬性的環(huán)境因子，雖然能通過環(huán)境因子擬合回歸方程對土壤屬性含量進行預測，但在鄉(xiāng)鎮(zhèn)尺度上環(huán)境因子的作用沒有縣域明顯，導致其逐步回歸方程的決定系數(shù)不高〔w(有機質(zhì))，R2=0.15；w(速效鉀)，R2=0.16〕(見表4)，空間制圖效果容易出現(xiàn)平滑現(xiàn)象，難以體現(xiàn)土壤屬性的空間變異性，因而不適合環(huán)境變化較單一的鄉(xiāng)鎮(zhèn)尺度(尤其平原地區(qū)). 而RK2法綜合了OK法和RK1法的優(yōu)點，通過把鄰近信息和環(huán)境因子同時納入逐步回歸方程，不管在預測精度(見表6)還是制圖效果上(見圖3)都優(yōu)于其他兩種方法，克服了OK法沒有考慮環(huán)境因子等宏觀因素和RK1法沒有考慮鄰近信息等局部因素的缺陷；RK2法較其他兩種方法在預測精度和揭示局部土壤屬性空間變異方面體現(xiàn)了一定的優(yōu)越性，空間分布模擬圖更加符合土壤屬性的空間變異性.

在鄉(xiāng)鎮(zhèn)尺度上，影響土壤屬性空間分布除定量因素(環(huán)境因子和鄰近信息)外，其他定性因素(如人為活動、土地利用類型、景觀類型、成土母質(zhì)、土屬、剖面構(gòu)型等)的影響也較為突出，引入更多的因素還可以進一步提高預測精度. 同時，多元回歸模型過于簡單，難以描述土壤屬性與自變量之間的復雜非線性映射關(guān)系，后續(xù)研究可應用具有較強的非線性逼近的神經(jīng)網(wǎng)絡來建立模型進行預測；也可以用其他插值方法(如集成神經(jīng)網(wǎng)絡、光滑薄板樣條等)來代替半方差函數(shù)進行殘差項OK法插值，使其即使在樣本不符合正態(tài)分布甚至極少樣本時也可以進行預測.

3 結(jié)論

a) 齊埠鎮(zhèn)耕地表層土壤w(有機質(zhì))、w(速效鉀)的平均值分別為35.03 g/kg和96.73 mg/kg，二者的自相關(guān)性較強，均為中等空間變異性. 相關(guān)性分析結(jié)果顯示，地形濕度指數(shù)和鄰近信息與兩種土壤屬性均有顯著相關(guān)，并且鄰近信息(P1)的相關(guān)系數(shù)達到0.6以上.

b) 通過引入環(huán)境因子和鄰近信息構(gòu)建RK2法、RK1法和OK法對齊埠鎮(zhèn)耕地表層土壤w(有機質(zhì))、w(速效鉀)空間分布進行模擬. 對62個采樣點進行建模，16個測試點進行獨立驗證的誤差分析表明，RK2法對土壤w(有機質(zhì))、w(速效鉀)預測結(jié)果的RMSE、MAE、MRE較OK法分別降低了18.06%、18.01%、21.77%和7.25%、9.49%、9.84%；較RK1法分別降低了22.48%、20.91%、22.02%和9.27%、12.61%、13.52%. 可以得出，僅基于環(huán)境因子的RK1法預測精度最低，基于環(huán)境因子和鄰近信息的RK2法預測效果最佳.

c) 由于該研究尺度為較小的鄉(xiāng)鎮(zhèn)尺度，環(huán)境因子較為單一，在一定程度上限制了僅基于環(huán)境因子的RK1法預測效果. 而引入鄰近信息后能夠明顯的提高土壤屬性的預測精度和空間分布模擬效果，能夠更加真實的反映和揭示土壤屬性空間變異規(guī)律. 但由于該研究選擇只考慮的環(huán)境因子和鄰近信息的影響，存在著一定的局限性，在后續(xù)研究中應充分考慮定量因素和定性因素還能進一步提高預測精度，進而為鄉(xiāng)鎮(zhèn)尺度上耕地利用和土壤環(huán)境保護提供更加準確的理論依據(jù).

[1] MCBRATNEY A B,PRINGLE M J.Estimating average and proportional variograms of soil properties and their potential use in precision agriculture[J].Precision Agriculture,1999,1(2):125- 152.

[2] DAI Wanhong,HUANG Yao.Relation of soil organic matter concentration to climate and altitude in zonal soils of China[J].Catena,2006,65(1):87- 94.

[3] SARKAR S,ROY A K,MARTHA T R.Soil depth estimation through soil-landscape modelling using regression Kriging in a Himalayan terrain[J].International Journal of Geographical Information Science,2013,27(12):2436- 2454.

[4] MEHNATKESH A,AYOUBI S,JALALIAN A,etal.Relationships between soil depth and terrain attributes in a semi arid hilly region in western Iran[J].Journal of Mountain Science,2013,10(1):163- 172.

[5] 劉二永,劉健,余坤勇,等.基于環(huán)境因子和R-STPS的林地土壤有機質(zhì)預測模型[J].農(nóng)業(yè)機械學報,2015,46(1):133- 137. LIU Eryong,LIU Jian,YU Kunyong,etal.Spatial prediction of forest soil organic matter based on environmental factors and R-ATPS interpolation methods[J].Transactions of the Chinese Society for Agricultural Machinery,2015,46(1):133- 137.

[6] KUMAR S,LAL R,LIU Desheng,etal.Estimating the spatial distribution of organic carbon density for the soils of Ohio,USA[J].Journal of Geographical Sciences,2013,23(2):280- 296.

[7] 馬泉來,高鳳杰,張志民,等.我國東北黑土丘陵區(qū)小流域土壤有機質(zhì)空間分布模擬[J].環(huán)境科學研究,2016,29(3):382- 390. MA Quanlai,GAO Fengjie,ZHANG Zhimin,etal.Simulation of spatial distribution of soil organic matter in a mollisol watershed in northeastern China[J].Research of Environmental Sciences,2016,29(3):382- 390.

[8] BHUNIA G S,SHIT P K,MAITI R.Spatial variability of soil organic carbon under different land use using radial basis function(RBF)[J].Modeling Earth Systems and Environment,2016.doi:10∶1007/s40808- 01s- 0070-x.

[9] 蘇曉燕,趙永存,楊浩,等.不同采樣點數(shù)量下土壤有機質(zhì)含量空間預測方法對比[J].地學前緣,2011,18(6):34- 40. SU Xiaoyan,ZHAO Yongcun,YANG Hao,etal.A comparison of predictive methods for mapping the spatial distribution of soil organic matter content with different sampling densities[J].Earth Science Frontiers,2011,18(6):34- 40.

[10] 史舟,李艷.地統(tǒng)計學在土壤學中的應用[M].北京:中國農(nóng)業(yè)出版社,2006.

[11] WANG Ku,ZHANG Chuanrong,LI Weidong.Comparison of geographically weighted regression and Regression Kriging for estimating the spatial distribution of soil organic matter[J].Giscience & Remote Sensing,2012,49(6):915- 932.

[12] ELDEIRY A A,GARCIA L A.Comparison of Ordinary Kriging,Regression Kriging,and Cokriging techniques to estimate soil salinity using LANDSAT images[J].Journal of Irrigation & Drainage Engineering,2010,136(6):355- 364.

[13] 張素梅,王宗明,張柏,等.利用地形和遙感數(shù)據(jù)預測土壤養(yǎng)分空間分布[J].農(nóng)業(yè)工程學報,2010,26(5):188- 194. ZHANG Sumei,WANG Zongming,ZHANG Bai,etal.Prediction of spatial distribution of soil nutrients using terrain attributes and remote sensing data[J].Transactions of the CSAE,2010,26(5):189- 194.

[14] WANG Binwu,ZHOU Weijun,MA Su,etal.Regression-kriging of soil organic matter using the environmental variables derived from MODIS and DEM[J].Agricultural Science & Technology,2012,13(4):838- 842.

[15] 李啟權(quán),王昌全,岳天祥,等.基于定性和定量輔助變量的土壤有機質(zhì)空間分布預測:以四川三臺縣為例[J].地理科學進展,2014,33(2):259- 269. LI Qiquan,WANG Changquan,YUE Tianxiang,etal.Prediction of distribution of soil organic matter based on qualitative and quantitative auxiliary variables:a case study in Santai County in Sichuan Province[J].Progress in Geography,2014,33(2):259- 269．

[16] LV Jianshu,LIU Yang,ZHANG Zulu,etal.Factorial Kriging and stepwise regression approach to identify environmental factors influencing spatial multi-scale variability of heavy metals in soils[J].Journal of Hazardous Materials,2013,261(13):387- 397.

[17] ZHAO Yongcun,XU Xianghua,DARILEK J L,etal.Spatial variability assessment of soil nutrients in an intense agricultural area,a case study of Rugao County in Yangtze River Delta Region,China[J].Environmental Geology,2009,57(5):1089- 1102.

[18] WERE K,SINGH B R,DICK? B.Spatially distributed modelling and mapping of soil organic carbon and total nitrogen stocks in the Eastern Mau Forest Reserve,Kenya[J].Journal of Geographical Sciences,2016,26(1):102- 124.

[19] 徐劍波,宋立生,夏振,等.基于GARBF神經(jīng)網(wǎng)絡的耕地土壤有效磷空間變異分析[J].農(nóng)業(yè)工程學報,2012,28(16):158- 165. XU Jianbo,SONG Lisheng,XIA Zhen,etal.Spatial variability of available phosphorus for cultivated soil based on GARBF neural network[J].Transactions of the CSAE,2012,28(16):158- 165.

[20] 李啟權(quán),王昌全,岳天祥,等.不同輸入方式下RBF神經(jīng)網(wǎng)絡對土壤性質(zhì)空間插值的誤差分析[J].土壤學報,2008,45(2):360- 365. LI Qiquan,WANG Changquan,YUE Tianxiang,etal.Error analysis of soil property spatial interpolation with RBF artificial neural network with different input methods[J].Acta Pedologica Sinca,2008,45(2):360- 365.

[21] 楊順華,張海濤,郭龍,等.基于回歸和地理加權(quán)回歸Kriging的土壤有機質(zhì)空間插值[J].應用生態(tài)學報,2015,26(6):1649- 1656. YANG Shunhua,ZHANG Haitao,GUO Long,etal.Spatial interpolation of soil organic matter using regression Kriging and geographically weighted regression Kriging[J].Chinese Journal of Applied Ecology,2015,26(6):1649- 1656.

[22] CHEN Chuanfa,ZHAO Na,YUE Tianxiang,etal.A generalization of inverse distance weighting method via kernel regression and its application to surface modeling[J].Arabian Journal of Geosciences,2014,8(9):1- 11.

[23] LIU Zhonggang,LI Manchun,SUN Yan,etal.Study on spatial autocorrelation of urban land price distribution in Changzhou City of Jiangsu Province[J].Chinese Geographical Science,2006,16(2):160- 164.

[24] 張甘霖,龔子同.土壤調(diào)查實驗室分析方法[M].北京:科學出版社,2012.

[25] 夏春蘭.多尺度條件下川中丘陵區(qū)土壤氮素分布的模擬研究[D].雅安:四川農(nóng)業(yè)大學,2013.

[26] HENGL T,GERARD B M,STEIN A,etal.A generic framework for spatial prediction of soil variables based on regression Kriging[J].Geoderma,2004,120(1/2):75- 93.

[27] CARLSON T N,RIPLEY D A.On the relation between NDVI,fractional vegetation cover,and leaf area index[J].Remote Sensing of Environment,1997,62(3):241- 252.

[28] 黃安,楊聯(lián)安,杜挺,等.基于多元成土因素的土壤有機質(zhì)空間分布分析[J].干旱區(qū)地理,2015,38(5):994- 1003. HUANG An,YANG Lianan,DU Ting,etal.Spatial distribution of the soil organic matter based on multiple soil factors[J].Arid Zone Research,2015,38(5):994- 1003.

[29] 張志霞,許明祥,師晨迪,等.黃土丘陵區(qū)不同地貌單元土壤有機碳空間變異的尺度效應[J].自然資源學報,2014,29(7):1173- 1184. ZHANG Zhixia,XU Mingxiang,SHI Chendi,etal.Scaling effect on spatial variability of soil organic carbon in different geomorphic units on the Loess Hilly Region[J].Journal of Natural Resources,2014,29(7):1173- 1184.

[30] 趙倩倩,趙庚星,姜懷龍,等.縣域土壤養(yǎng)分空間變異特征及合理采樣數(shù)研究[J].自然資源學報,2012,30(1):31- 35. ZHAO Qianqian,ZHAO Gengxing,JIANG Huailong,etal.Study on spatial variability of soil nutrients and reasonable sampling number at county scale[J].Journal of Natural Resources,2012,30(1):31- 35.

[31] HERNANDEZSTEFANONI J L,PONCEHERNANDEZ R.Mapping the spatial variability of plant diversity in a tropical forest:comparison of spatial interpolation methods[J].Environmental Monitoring & Assessment,2006,117(1/2/3):307- 334.

[32] ZHU Qu,LIN Haisheng.Comparing Ordinary Kriging and Regression Kriging for soil properties in contrasting landscapes[J].Pedosphere,2010,20(5):594- 606.

[33] 王彬武,周衛(wèi)軍,馬蘇,等.耕地土壤速效鉀含量的空間預測方法研究[J].地理與地理信息科學,2011,27(5):91- 94. WANG Binwu,ZHOU Weijun,MA Su,etal.Spatial prediction method of available potassium contents in arable soil[J].Geography and Geo-Information Science,2011,27(5):91- 94.

Prediction of Spatial Distribution of Soil Properties Based on Environmental Factors and Neighbor Information

JIANG Yefeng1,2, SUN Kai1,2, GUO Xi1,2*, YE Yingcong2, RAO Lei1,2, LI Weifeng1,2

1.Academy of Land Resource and Environment, Jiangxi Agricultural University, Nanchang 330045, China 2.Key Laboratory of Poyang Lake Watershed Agricultural Resources and Ecology of Jiangxi Province, Nanchang 330045, China

2016-12-27

2017-02-21

國家自然科學基金項目(41361049)；江西省自然科學基金項目(20122BAB204012)；江西省贛鄱英才“555”領(lǐng)軍人才項目(201295)

江葉楓(1994-), 男, 江西余干人, jiangyf0308@163.com.

X142

1001- 6929(2017)07- 1059- 10

A

10.13198/j.issn.1001- 6929.2017.02.09

*責任作者，郭熙(1974-), 男, 江西永豐人, 副教授, 博士, 主要從事土壤遙感與信息技術(shù)及土地資源利用研究, xig435@163.com

江葉楓,孫凱,郭熙,等.基于環(huán)境因子和鄰近信息的土壤屬性空間分布預測[J].環(huán)境科學研究,2017,30(7):1059- 1068.

JIANG Yefeng,SUN Kai,GUO Xi,etal.Prediction of spatial distribution of soil properties based on environmental factors and neighbor information[J].Research of Environmental Sciences,2017,30(7):1059- 1068.