基于協同克里格的耕層土壤速效鉀空間異質性研究

江葉楓，郭 熙

(江西農業大學 江西省鄱陽湖流域農業資源與生態重點實驗室，江西 南昌 330045)

土壤是有限且不可再生的自然資源，在提供生態系統服務方面發揮著至關重要的作用[1-2]。為維護自然資源以支持未來的可持續發展，亟須改善和維護土壤的生態系統服務功能[3]。土壤屬性的變化與土壤生態系統服務功能的可持續性息息相關，土壤侵蝕、土壤酸化、土壤重金屬污染與土壤養分供需不平衡等，都不可避免地會影響土壤質量和生態系統服務功能，降低土壤對人類的支持能力[3-5]。土壤養分是土壤可持續利用和管理的重要基礎[3-4]。土壤速效鉀(AK)能夠直觀反映土壤鉀素水平和供應能力，作為作物生長周期內可以獲取的主要養分，直接影響農產品產量和品質[6]。近年來，由于各地區生產水平、耕作技術措施、種植制度和耕地資源利用強度不同，土壤速效鉀的空間異質性普遍存在[7]。隨著精準農業和配方施肥技術的發展，土壤速效鉀的空間異質性研究不斷深入[8]。準確掌握土壤速效鉀的空間異質性對科學制定配方施肥方案、保護土壤、實現土壤的可持續利用等具有重要意義[6-8]。

土壤是連續的，科學家難以檢測到每一處土壤的屬性值，不得不依賴于空間插值技術[9]。國內外學者主要運用克里格空間插值法對土壤屬性進行空間插值并研究土壤屬性的異質性[10-11]。該方法基于半變異函數，是對土壤屬性進行的最優無偏估測，能有效揭示土壤屬性空間分布和變異特征。其中，普通克里格法(ordinary Kriging，OK)只需要考慮目標變量，簡單易行，應用廣泛，但插值結果在土壤屬性異質性描述方面受采樣數量和密度的影響[7，12]。協同克里格法(co-Kriging，COK)是一種基于目標變量與輔助變量的半變異函數，由于輔助變量與土壤屬性相關性較高，因而以這種方式生成的可視化圖通常能夠更加高效、客觀、真實地描述土壤屬性空間變異特征，在土壤屬性空間變異研究中越發受到研究人員的青睞[13-14]。

COK較OK能提高空間插值精度，但在土壤屬性的異質性描述方面是否也較OK存在一定優勢尚不清楚[13-14]。當下，大部分研究對土壤屬性的異質性描述仍處于定性分析層面，缺乏對影響因子作用的量化分析[7]。從相關文獻來看，關于土壤性質變異的影響因素定量分析報道仍存在明顯不足：1)對土壤屬性的輔助變量考慮不充分，地形因子因能較好地刻畫土壤性質的地帶性分異規律而受到較大關注[12，15]，但僅少量學者定量分析了成土母質和土壤類型等對土壤屬性空間異質性的影響[11，16]；2)相關的定量研究較多針對土壤多種性質或針對有機碳[17]、氮素[4]和磷素[18]等單一性質，鮮見關于土壤鉀素的報道，特別是在我國耕地復種指數增高和土壤鉀素供應嚴峻的背景下，對速效鉀變異的影響因素缺乏量化認識[7]；3)不同區域COK輔助變量的選取各有不同，前人曾對川中丘陵[19]、黃土高原[20]、西北山區[21]和城鄉交錯帶[22]等進行研究，而鄱陽湖平原區則鮮見報道。

為此，本研究以鄱陽湖平原區典型縣——萬年縣為案例區，運用測土配方施肥采集的耕層(0～20 cm)土壤樣品數據，以土壤速效鉀為研究對象，量化土壤速效鉀空間變異的影響因素，運用OK和COK對土壤速效鉀的空間異質性進行研究。結果可為類似地區土壤鉀素調控、土地可持續管理和利用，以及土壤屬性空間分布預測模型構建提供參考。

1 材料與方法

1.1 研究區概況

研究區位于萬年縣，屬鄱陽湖平原區，地理坐標介于116°46′—117°15′E、28°30′—28°54′N之間，總面積1 140.76 km2，東西長47 km，南北寬43 km。研究區屬亞熱帶季風性氣候，地貌類型以平原為主，呈“西北部高、東南部低”的總體趨勢(圖1-a)，海拔12～651 m。土類主要是紅壤和水稻土，并伴隨有極少量的潮土、石灰土和紫色土，亞類主要為紅壤和潴育型水稻土，土屬主要包括黃泥紅壤和鱔泥紅壤，潴育型潮砂泥田、潴育型黃泥田和潴育型鱔泥田。成土母質包括紫色泥頁巖類風化物、第四紀紅色黏土、河流沖積物、泥質巖類風化物、酸性結晶巖類風化物、碳酸巖類風化物和紫色砂礫巖類風化物。境內水資源比較豐富，河流主要包括樂安河、珠溪河和萬年河等182條，總長806 km，河網密度707 m·km-2，年均降水量1 766 mm，年均無霜期263 d，年均相對濕度82%。土地利用分布格局為“六山一水二分田，一分道路和莊園”，耕地面積約為230 km2，耕層土壤pH介于4.5～7.0。萬年縣是世界稻作文化的發源地，也是人工栽培稻起源地和貢米的原產地，被稱為“稻米之鄉”，傳統的稻米習俗已在該地傳承了幾千年，萬年貢米是晚秈稻“塢源早”品種加工而成的產品，原只產于萬年縣裴梅鎮部分鄉村，后經種植推廣，截至2016年，萬年縣貢米種植面積達170 km2。

1.2 土壤采樣與數據處理

于2014年10—11月在作物收割后，根據《測土配方施肥技術模式》，在充分考慮土壤類型、成土母質和地形條件的基礎上，遵循均勻性、代表性和連續性的原則，運用“S”形采樣策略和多點混合的方法采集耕層(0～20 cm)土壤樣品，每個樣點采集1 kg土樣，4個重復樣本在該點周圍5 m范圍內采集，并運用GPS詳細記錄該點的經緯度、成土母質、土壤類型和海拔等信息，共獲取100個耕層土壤樣品。土壤樣品經自然風干后，帶回實驗室磨碎過篩，采用乙酸銨浸提—原子吸收分光光度法測定土壤速效鉀含量[7]。數字高程模型(digital elevation model， DEM)如圖1-a所示，分辨率30 m，通過地理空間數據云(http://www.gscloud.cn/)獲取。

受采樣及化學分析誤差的影響，土壤速效鉀含量測量結果存在離群值。采用拉依達準則法[4]對測定數據進行檢驗，剔除離群值，得到剔除后有效樣點96個(圖1-b)。高程、坡度、坡向、曲率、坡度變率和坡向變率均由數字高程模型通過ArcGIS 10.2表面分析模塊進行分析和提取，各地形因子的計算公式詳見文獻[23]。土壤AK根據第二次土壤普查土壤養分分析標準分為6級，1～6級分別對應極高(>200 mg·kg-1)、豐富(>150～200 mg·kg-1)、中等(>100～150 mg·kg-1)、較低(>50～100 mg·kg-1)、低(>30～50 mg·kg-1)和極低(≤30 mg·kg-1)。成土母質和土壤類型為定性變量，采用啞變量[4，16]進行賦值。

1.3 模擬精度評價

采用交叉驗證評價OK和COK的預測精度。交叉驗證利用每個實測點周圍點對該點進行預測，將該點預測值與實測值進行比較。以均方根誤差(RMSE)、平均絕對誤差(MAE)、平均相對誤差(MRE)對預測值與實際采樣值進行對比分析，得出精度評價結果，計算公式如下：

(1)

(2)

(3)

圖1 研究區數字高程模型和樣點分布圖Fig.1 Maps of DEM and sample points in study area

描述性統計分析、相關性分析和回歸分析在IBM SPSS statistics 22.0中實現，半變異函數分析在GS+9.0中實現，OK和COK插值通過ArcGIS 10.2實現。

2 結果與分析

2.1 描述性統計分析

2.1.1 土壤速效鉀

研究區土壤速效鉀值域范圍為33.46～164.84 mg·kg-1，均值為82.04 mg·kg-1。根據第二次土壤普查土壤養分分級標準，研究區速效鉀含量為4級，屬于較低水平，勉強能滿足當地耕地生產需要。因部分地區土壤鉀素盈余或不足，以及耕作技術措施差異，同時考慮耕地復種指數和化肥施用比例不協調等原因，建議重視土壤鉀素保育[7]。研究區土壤速效鉀的最大值與最小值之比為4.93，表明研究區土壤鉀素分布極不平衡，變異系數為30.49%，屬于中等變異。Kolmogorov-Smirnov(K-S)檢驗的結果表明，土壤速效鉀含量符合正態分布。

2.1.2 土壤類型與成土母質

土壤類型反映區域條件的綜合變化，因成土過程、礦物組成和發育程度不同，土壤速效鉀存在空間異質性[7，16]。由表1可知，不同類型土壤速效鉀均值差異顯著(P<0.05)，表現出水稻土>紅壤的總體趨勢。這主要是由于水稻土受長期人為活動施鉀肥的影響[4]，導致其速效鉀均值要高于紅壤。各土屬中，以潴育型潮砂泥田均值最高，達到100.88 mg·kg-1，鱔泥紅壤最低，為73.95 mg·kg-1。

成土母質的礦物分解是土壤鉀素的基本來源，不同成土母質因土壤團聚體數量及其穩定性、物理化學組成和風化淋溶進程等引起土壤速效鉀變異[24]。從表2可以看出，成土母質對土壤速效鉀有顯著(P<0.05)影響。不同母質間土壤速效鉀含量為61.10～92.04 mg·kg-1，其中，以白堊紀紫色泥頁巖類風化物發育而來的土壤速效鉀均值含量最高，而紫色砂礫巖類風化物最低。白堊紀紫色泥頁巖一般含碳酸鈣，呈中性或微堿性反應，有機質含量低，但磷、鉀豐富；酸性結晶巖類風化物以殘積物和坡積物為主，質地主要為黏壤土，土壤保水保肥能力強，土壤速效鉀含量也相對較高；紫色砂礫巖類風化物易風化，水土流失嚴重，因此土壤速效鉀較低[7，24]。其他研究也得到類似結論[16]。

表1不同土壤類型速效鉀描述性統計特征

Table1Descriptive statistic characteristics of soil available potassium relative to soil types

土類Soil group亞類Subgroup土屬Soil family樣點數Sample No.均值Mean/(mg·kg-1)SD/(mg·kg-1)CV/%紅壤Red soil紅壤Red soil黃泥紅壤Yellow mud red soil584.84 b27.0131.84鱔泥紅壤Muddy red soil3673.95 c18.3624.83水稻土潴育型水稻土潴育型潮砂泥田Breeding type tidal mud field8100.88 a29.7629.50Paddy soilWaterloggogenic潴育型黃泥田Breeding yellow mud field1489.61 b20.8023.21type paddy soil潴育型鱔泥田Breeding muddy field3382.66 b27.5233.29

SD，標準差；CV，變異系數。同列數據后無相同字母的表示差異顯著(P<0.05)。下同。

SD， Standard deviation; CV， Coefficient variation. Data marked without the same letters indicated significant difference atP<0.05. The same as below.

表2不同成土母質土壤速效鉀描述性統計特征

Table2Descriptive statistic characteristics of soil available potassium relative to parent materials

成土母質Parent material樣點數Sample No.均值Mean/(mg·kg-1)SD/(mg·kg-1)CV/%紫色泥頁巖類風化物Purple mud shale weathering1492.04 a29.6232.18第四紀紅色黏土Quaternary red clay2989.72 a26.2029.20河流沖積物River alluvial1081.78 b14.6417.90泥質巖類風化物Mudstone weathering1269.57 b13.6919.68酸性結晶巖類風化物Acidic crystalline rock weathering589.74 a19.9022.18碳酸巖類風化物Carbonate weathering2175.15 b24.3232.36紫色砂礫巖類風化物Purple glutenite weathering561.10 c10.8917.82

2.1.3 土壤速效鉀與影響因素的相關性分析

從表3可以看出，土壤速效鉀與pH的相關性達極顯著水平(P<0.01)，相關系數為-0.258，說明在一定范圍內土壤速效鉀隨pH降低而逐漸增加，這與其他區域研究結果一致[25]。土壤陽離子交換量(CEC)與速效鉀相關性同樣達極顯著水平(P<0.01)，且相關系數達到0.602，表明CEC與AK呈較強的正相關關系，CEC越高，土壤速效鉀越容易累積。高程與土壤速效鉀的相關性達極顯著水平(P<0.01)，相關系數為-0.414，高程越高處土壤受到的沖刷侵蝕越嚴重，土壤速效鉀越易流失，而地勢低洼處速效鉀易隨地表物質積聚[7]。土壤速效鉀與其他地形因子，如坡度、坡向、曲率、坡度變率、坡向變率的相關性不顯著，可能與鄱陽湖平原區有關，地形因子本身的變異較小。同時，本研究獲取的DEM數據精度較低，其派生的地形變量精度也相對較低，這也會在一定程度上影響地形因子與土壤速效鉀的相關性[16]。

表3土壤速效鉀與pH、CEC和地形因子的Pearson相關系數

Table3Pearson correlation coefficients of soil available potassium with pH， CEC， and terrain factors

指標IndexrpH-0.258**CEC0.602**高程Elevation-0.414**坡度Slope-0.098坡向Aspect-0.072曲率Curvature0.18坡度變率Slope variability-0.073坡向變率Aspect variability-0.06

*，P<0.05; **，P<0.01.

2.1.4 土壤速效鉀與影響因素的回歸分析

為定量揭示成土母質、土壤類型(亞類和土屬)、高程、pH和CEC對研究區耕層土壤速效鉀空間變異的影響，對影響因素進行回歸分析(表4)。結果表明：成土母質、土壤類型(亞類和土屬)、高程、pH和CEC對土壤速效鉀空間變異影響均達極顯著水平(P<0.01)。成土母質對土壤速效鉀空間變異的影響程度較低，為8.3%，這與相關研究結果[16]一致。江葉楓等[12，16]研究表明，該區主要成土母質包括白堊紀紫色泥頁巖、第四紀紅色黏土和幾種風化物，多由花崗巖、石灰巖和紫色砂巖經風化而成，巖性總體上較為一致，導致土壤機械組成相近[16]，成土母質主要通過影響土壤機械組成來影響土壤速效鉀分異，因此，成土母質對本研究區土壤速效鉀影響程度較低。土屬的影響程度要高于土類和亞類。江葉楓等[26]認為，由于土屬較土類和亞類而言反映的成土過程、化學成分和生物活性等信息更多，因此影響程度更高。研究區pH介于4.5～7.0，當pH值降低時，土壤膠體微粒表面所負電荷也減少，導致K+吸附量隨之增加。陳洋等[7]研究表明，在酸性環境中K+運移量較少，固定量增多。CEC基本上代表了土壤可能保持的養分數量，其影響因素包括土壤膠體類型、土壤質地和pH值，這些都是影響土壤速效鉀變異的因素。因此，CEC對土壤速效鉀變異影響程度較高。

2.2 土壤速效鉀的半變異函數分析

表5給出了土壤速效鉀半變異函數擬合結果。

表4不同因素對土壤速效鉀的回歸分析

Table4Regression analysis of soil available potassium with affecting factors

影響因素FactorF決定系數R2調整決定系數Adjusted R2P成土母質 Parent material1.9130.0920.083<0.01土壤類型 Soil group土類 Soil group5.4180.0710.066<0.01亞類 Subgroup5.4180.0710.066<0.01土屬 Soil family6.2980.0850.080<0.01高程 Elevation12.6900.1690.165<0.01pH4.9060.0650.059<0.01CEC30.3760.3560.354<0.01

表5土壤速效鉀的半變異函數值

Table5Semivariance parameters for soil available potassium

土壤屬性Soil properties模型Model塊金值Nugget基臺值Partial sill塊金效應Ratio of nuggest to sill/%變程Range/km決定系數Determination coefficient殘差Residual速效鉀AK指數Exponential260.54560.0046.532.640.6581.37×10-4

基于最大的擬合系數和最小的殘差得到土壤速效鉀的最佳擬合模型為指數模型，塊金值為260.54，說明存在由實驗誤差和田間采樣等人為因素造成的空間變異，塊金效應值為46.53%，表明速效鉀的空間變異受自然特征和人為活動的共同影響，為中等程度的空間相關性[16-20]。

2.3 土壤速效鉀空間分布模擬

為直觀反映研究區土壤速效鉀的空間異質性，在半變異函數的基礎上運用OK、COK1、COK2和COK3對研究區土壤速效鉀進行空間插值，用可視化圖作為手段來評估4種方法描述土壤速效鉀空間異質性的能力。其中，COK1代表以CEC作為輔助變量進行插值，COK2代表以CEC和pH作為輔助變量進行插值，COK3代表以CEC、pH和高程同時作為輔助變量進行插值。

交叉驗證的精度評價結果(表6)表明，基于輔助變量的3種COK明顯優于僅基于目標變量的OK，COK1、COK2和COK3較OK的RMSE分別降低了1.03、1.92、4.86 mg·kg-1，MAE分別降低了1.10、2.21、5.37 mg·kg-1，MRE分別降低了1.41、2.74、5.50個百分點。基于輔助變量的COK預測精度得到了較為明顯的提高。從表6還可以看出，隨著輔助變量增加，COK的精度也在相應提升。這一方面表明增加與目標變量相關的輔助變量協助空間插值可以提高精度，另一方面也說明運用多個輔助變量較單個輔助變量精度更高，后期研究中可以適當增加輔助變量個數以提高預測精度。

從圖2可以看出，OK、COK1、COK2和COK3對研究區土壤速效鉀的模擬均表現出“西北部高、東南部低”的總體趨勢，這與數字高程模型的空間變化較吻合，比較符合地學分布規律。從空間分布模擬效果看，4種方法模擬的局部特征差異明顯。OK得到的空間分布模擬結果較平滑，高低值界限較清晰，難以準確地表達土壤速效鉀的空間異質性，預測值域范圍58.83～121.83 mg·kg-1，與統計分析值有較大差距。COK1運用相關性最強的CEC作為輔助變量進行空間插值，預測精度較OK有所提升，但空間分布模擬效果提升不明顯，從預測范圍(63.59～109.48 mg·kg-1)和均值(78.59 mg·kg-1)來看，COK1預測更趨向于均值，平滑效應明顯，難以較清晰地刻畫土壤速效鉀的空間異質性信息。COK2在COK1的基礎上引入pH作為輔助變量進行協同插值，不僅預測精度有所提升，在空間異質性信息描述方面也有很大提升。COK2預測土壤速效鉀值域介于58.19～122.08 mg·kg-1，比較接近統計分析值，預測的空間分布模擬圖高低值呈塊狀分布，同時出現了較多高值區域包含低值或低值區域包含高值部分，能更詳細地描述土壤速效鉀的空間異質性。COK3預測結果高低值呈塊狀分布，空間連續性較其他3種方法有所增強，插值精度明顯提高，能刻畫更多速效鉀空間異質性的細節信息，土壤速效鉀為突變而非漸變，且預測范圍為36.78～124.31 mg·kg-1，與描述性數據分析結果最為接近。

表6不同方法精度對比

Table6Precision comparison of different methods

方法MethodsRMSE/(mg·kg-1)MAE/(mg·kg-1)MRE/%OK25.1119.2625.23COK124.0818.1623.82COK223.1917.0522.49COK320.2513.8919.73

3 討論

本研究表明，研究區土壤速效鉀普遍較低，空間分布受自然特征和人為活動及其協同作用的共同影響，但空間變異主要是由人為活動導致的。Pearson相關性分析結果表明pH、CEC和高程與土壤速效鉀相關性達極顯著水平。單因素方差分析結果表明，不同成土母質和土壤類型的土壤速效鉀差異顯著，不同因素的影響程度由大到小依次為CEC>高程>成土母質>土壤類型(亞類和土屬)>pH。以3個輔助變量協同進行空間插值的COK3模擬精度最高，空間分布表現出“西北部高、東南部低”的總體趨勢，空間異質性描述方面更加符合研究區實際情況和地學分布規律。

圖2 不同方法的土壤速效鉀空間分布預測結果Fig.2 Maps of soil available potassium by different methods

研究區耕層土壤速效鉀均值為82.04 mg kg-1，含量處于較低水平，與孫凱等[27]、江葉楓等[28]和劉雪梅等[29]對鄱陽湖平原區的研究結果較為接近。從不同地貌類型看，均值要低于川西山區耕地[21]、秦嶺山地區[30]和黃土高原區[20]。從不同土地利用類型來看，也要低于川中丘陵植煙區[25]、江南茶區[31]和冬小麥-夏玉米輪作區[32]，同時也低于福建省耕地[33]和江蘇省農田[34]。這表明研究區土壤速效鉀的分異可能是由于地貌類型和土地利用方式引起的。鄱陽湖平原區水熱資源豐富，基巖風化、分解速率快，淋溶作用強烈，同時作為中國的商品糧基地，其土地利用強度也必然較大，土壤中鉀素損失嚴重。

由于土壤屬性獲取的昂貴性與費時性，為了得到精確的土壤屬性空間分布，應充分考慮影響土壤屬性空間分布的因素及土壤屬性的空間自相關性[28]。目前，地統計學中的克里格插值方法基于半變異函數的結構性，在考慮空間自相關的基礎上對未知采樣區域的土壤屬性進行無偏最優估值，廣泛應用于土壤屬性的空間異質性研究[10-16]。OK只簡單地考慮被預測土壤屬性的空間自相關，僅關注于目標變量，而沒有充分利用影響土壤屬性空間分布的各種因素，相比之下，COK可以利用與目標變量相關性較好的輔助變量來提高目標變量的估測精度，能充分考慮影響目標變量的其他因子，將目標變量的空間自相關性和輔助變量間的交互相關性結合起來用于無偏最優估計，因而更有利于提高估測精度，這與前人研究結果一致[13-14]。本研究還發現，對于平原區，單個的輔助變量可能會在土壤屬性異質性描述方面存在一定的偏見，而增加輔助變量的個數可以更加詳細地刻畫其異質性，在后期研究中應注意這一現象。

一般認為，耕地土壤速效養分含量受農業生產活動影響較大[7]。半變異函數分析結果表明，研究區土壤速效鉀空間變異受自然特征和人為活動及其協同效應的共同作用[4]。成土母質和土壤類型是土壤速效鉀空間分異的自然特征，對土壤速效鉀空間變異影響程度顯著。黏土礦物是土壤鉀素供應的初始來源。形成各種土壤類型的基巖、黏土礦物、成土母質等類型不盡相同，母質遺傳了基巖礦物特性，并構成了鉀素有效性的載體[35]。盡管土壤鉀素的供應潛力受制于自然特征，但自然特征在一定時期內不會發生顯著變化[4]，因此研究區土壤速效鉀空間變異主要是由于人為活動導致的。施肥、灌溉、翻耕，以及復種指數等均會引起土壤速效鉀空間變異，然而輸入與輸出并非均衡，因此不同地區含量也有所差異。

本研究基于COK3統計表明，AK含量處于4級的面積為87.41%，而第二次土壤普查結果顯示3級及以上的面積為55.97%，表明研究區絕大部分耕地土壤速效鉀含量下降程度較大。考慮到鄱陽湖平原區是商品糧基地，為維持土壤和農業可持續發展，建議調控復種指數，增加有效鉀施用。