典型白漿土區(qū)耕層土壤速效鉀空間異質(zhì)性研究

侯永華，王 霞，趙映慧，羅 沖，劉煥軍，馬士耐，王 爍

（1. 東北農(nóng)業(yè)大學(xué) 公共管理與法學(xué)院，黑龍江 哈爾濱 150030；2. 中國科學(xué)院 東北地理與農(nóng)業(yè)生態(tài)研究所，吉林 長春 130012）

白漿土主要分布在黑龍江省和吉林省的東北部[1]，是兩省的主要耕作土壤類型之一。其中，黑龍江省白漿土面積約為331.2 hm2，主要分布于該省的三江平原地區(qū)。白漿土常與黑土、草甸土等混雜分布，三者中黑土層在相鄰?fù)寥乐凶顪\薄，厚度介于6～20 cm，平均為12.8 cm[2]。白漿土成土母質(zhì)主要是第四紀(jì)河湖沉積物，有一層瘠薄的白漿層，質(zhì)地黏重，且透水性差[3]。因此，該類型土壤表層旱澇頻發(fā)，且在長期耕作后易造成土壤偏酸性，對作物生長發(fā)育產(chǎn)生不良影響。土壤養(yǎng)分是土壤可持續(xù)利用與管理的重要基礎(chǔ)，為作物生長發(fā)育提供了各種營養(yǎng)元素，是不可再生資源[4-5]。土壤速效鉀（Available kalium，AK）作為作物生長過程中可獲取的主要養(yǎng)分之一，直接影響著農(nóng)產(chǎn)品的質(zhì)量與品質(zhì)[6]。預(yù)測白漿土區(qū)AK 含量及其空間分布，對于白漿土改良、糧食產(chǎn)量提升等具有重要意義。

國內(nèi)外學(xué)者主要利用克里格空間插值法對土壤屬性的異質(zhì)性進(jìn)行研究[7-8]。伴隨精準(zhǔn)農(nóng)業(yè)和測土配方技術(shù)的應(yīng)用與發(fā)展，學(xué)者們對土壤全氮（Total nitrogen，TN）[9-10]、有 效 磷（Available phosphorus，AP）[11-12]、有 機(jī) 質(zhì)（Organic matter，OM）[13-14]等養(yǎng)分及微量元素[15]開展了定量研究，以分析土壤的空間變異特性。克里格插值可實(shí)現(xiàn)未采樣區(qū)域土壤養(yǎng)分變化量的無偏最優(yōu)估計(jì)，常被應(yīng)用于土壤養(yǎng)分空間異質(zhì)性研究。普通克里格（Ordinary kriging，OK）是GIS（地理信息系統(tǒng)）地統(tǒng)計(jì)分析中重要的分析方法之一。與OK 相比，協(xié)同克里格（Cokriging，COK）可添加與目標(biāo)變量相關(guān)的其他環(huán)境輔助變量，根據(jù)與土壤養(yǎng)分相關(guān)性強(qiáng)弱實(shí)現(xiàn)諸多變量的疊加分析，削弱樣本采集過程中如環(huán)境、人工等其他因素的影響，從而更加客觀和全面地反映土壤養(yǎng)分的空間異質(zhì)性。因此，COK 成為當(dāng)前土壤養(yǎng)分研究常用的分析方法之一。石淑芹等[16]的研究表明，通過添加區(qū)域協(xié)同變量的COK 插值，可以提高土壤養(yǎng)分空間預(yù)測準(zhǔn)確度；黃安等[17]的研究表明，控制輔助變量進(jìn)入COK 插值順序有助于提升模型預(yù)測精度；章清等[18]基于主成分分析篩選了輔助變量，并以此提升了插值模型的擬合精度；馮家東等[19]認(rèn)為，篩選相關(guān)性高的輔助變量更有利于反映土壤養(yǎng)分的空間異質(zhì)性。所研究區(qū)域不同，往往選取的輔助變量會(huì)有不同，但總體而言，COK 插值的擬合效果明顯好于OK 插值法。目前尚未有白漿土區(qū)域速效鉀空間異質(zhì)性相關(guān)研究，因此，擬以黑龍江省樺南縣曙光農(nóng)場（46°13′～48°23′N，130°17′～130°39′E；典型白漿土區(qū)域）為研究區(qū)，通過分析和量化AK 空間影響因素，并將與AK 相關(guān)性較強(qiáng)的輔助因子作為輔助環(huán)境變量進(jìn)行COK 插值，研究典型白漿土區(qū)耕層土壤AK 的空間變異規(guī)律，為典型白漿土區(qū)土壤改良提供理論依據(jù)，進(jìn)一步提升典型白漿土區(qū)作物產(chǎn)量，保障糧食安全。

1 材料和方法

1.1 研究區(qū)概況

曙光農(nóng)場位于佳木斯市樺南縣，屬于三江平原地區(qū)。場區(qū)東西長35 km，土地總面積177 km2，其中，耕地面積超1萬hm2。地勢北高南低，東高西低，呈現(xiàn)為由東向西傾斜的魚脊形。土壤類型以白漿土為主，并有少部分黑土、草甸土分布，成土母質(zhì)主要是第四紀(jì)河湖沉積物。該農(nóng)場主要種植作物有玉米、大豆、水稻等。

1.2 土壤采樣與土壤數(shù)據(jù)處理

綜合曙光農(nóng)場作物類型、土壤類型和地形等不同因素，在遵循全面性、均衡性和代表性的基礎(chǔ)上布設(shè)采樣點(diǎn)，利用S 形取樣法采集0～20 cm 耕層土壤。采樣時(shí)間為2019 年10 月，共采集土樣186 份（圖1）。在采樣時(shí)記錄采樣點(diǎn)經(jīng)緯度、高程（Elevation，E）等信息，帶回實(shí)驗(yàn)室晾干，研磨，過0.149 mm 篩后測pH 值、OM 含量、AP 含量等。AK含量用乙酸銨浸提-火焰光度法測定；pH 值用電位法測定；OM 含量采用油浴加熱重鉻酸鉀氧化-容量法測定；AP 含量采用Olsen 法測定;AN（堿解氮）含量使用堿解擴(kuò)散法測定。該農(nóng)場DEM（數(shù)字高程模型，圖2）數(shù)據(jù)分辨率為30 m，數(shù)據(jù)來源于美國地質(zhì)勘探局網(wǎng)站（http：//www.usgs.gov/），坡度、坡向等信息利用DEM數(shù)據(jù)通過ArcGIS 10.2軟件提取。

圖1 典型白漿土區(qū)耕層土壤AK實(shí)地采樣點(diǎn)分布Fig.1 AK field sampling point distribution of topsoil in typical albic soil area

圖2 曙光農(nóng)場DEMFig.2 Shuguang farm DEM

土壤數(shù)據(jù)通過SPSS 26.0 進(jìn)行描述性統(tǒng)計(jì)分析和正態(tài)分布檢驗(yàn)（K-S 檢驗(yàn)），利用GS 9.0+軟件對AK 進(jìn)行半變異函數(shù)分析和理論模型擬合，利用ArcGIS 10.2軟件中地統(tǒng)計(jì)模塊OK插值和COK插值進(jìn)行空間分布特征分析。主要使用皮爾森相關(guān)系數(shù)對AK 的相關(guān)因素進(jìn)行分析，在此基礎(chǔ)上選擇最佳輔助變量用于COK插值。

1.3 地統(tǒng)計(jì)分析

地統(tǒng)計(jì)學(xué)利用原始數(shù)據(jù)和半方差函數(shù)的結(jié)構(gòu)性對未采樣點(diǎn)進(jìn)行無偏差估計(jì)。本研究主要采用地統(tǒng)計(jì)學(xué)中的半變異和克里格插值法研究AK 的空間變異性，主要公式：

式中，γ(h)為間距的半方差，在一定范圍內(nèi)隨x（兩樣本間的距離）的增加而增大，當(dāng)測點(diǎn)間距大于最大相關(guān)距離時(shí)，該值趨于穩(wěn)定。Q(x)、Q(x+h)分別為隨機(jī)變量Q在空間位置x、x+h上的取值。

1.4 遙感數(shù)據(jù)獲取與處理

獲取裸土?xí)r期（2019 年10 月3 日）Sentinel-2A影 像（https：//code.earthengine.google.com/）數(shù) 據(jù)。Sentinel-2A 影像為幾何校正和大氣處理后的地表反射率數(shù)據(jù)，Sentinel-2A 通過SNAP（哨兵應(yīng)用平臺）重采樣為10 m 空間分辨率。在此基礎(chǔ)上計(jì)算該時(shí)期歸一化差異植被指數(shù)（Normalized difference vegetation index，NDVI）、歸 一 化 差 分 水 體 指 數(shù)（Normalized difference water index，NDWI）、比值植被指數(shù)（Ratio vegetation index，RVI）等。將這些遙感指數(shù)與AK 含量進(jìn)行相關(guān)性分析后擇優(yōu)選取作為輔助變量進(jìn)行COK 插值，將所得結(jié)果與OK 插值所得結(jié)果及其他COK插值所得結(jié)果進(jìn)行對比分析。

1.5 擬合精度評價(jià)

通過交叉檢驗(yàn)對OK和COK插值所得結(jié)果進(jìn)行驗(yàn)證。主要選取均方根誤差（RMSE）、平均相對誤差（MRE）和平均絕對誤差（MAE）作為檢驗(yàn)精度的評價(jià)指標(biāo)，計(jì)算公式：

式中，X*i為樣點(diǎn)的預(yù)測值；Xi為樣點(diǎn)實(shí)測值；n為樣點(diǎn)數(shù)。RMSE、MRE和MAE值越小，表示插值擬合效果越好，預(yù)測精度越高。

2 結(jié)果與分析

2.1 典型白漿土區(qū)耕層土壤AK描述性統(tǒng)計(jì)分析

研究區(qū)AK 含量介于92.00～364.00 mg/kg，平均為171.78 mg/kg（表1）。為使該數(shù)據(jù)符合正態(tài)分布，將其進(jìn)行對數(shù)變換。變換后AK 含量介于1.96～2.56 mg/kg，平均為2.22 mg/kg。AK含量變異系數(shù)由27.52%降低為對數(shù)變換后的5.32%。根據(jù)全國第2次土壤普查標(biāo)準(zhǔn)[AK含量＞100.00 mg/kg處于三級高水平]，本研究區(qū)AK含量大部分處于高水平范圍內(nèi)。AK 含量最大值與最小值差異較大，表明AK 含量在研究區(qū)內(nèi)的值域分布較廣。

表1 典型白漿土區(qū)耕層土壤AK總體描述性統(tǒng)計(jì)Tab.1 Overall descriptive statistics of AK in typical albic soil area

2.2 典型白漿土區(qū)耕層土壤AK主要相關(guān)因素分析

2.2.1 典型白漿土區(qū)耕層土壤AK 輔助變量相關(guān)性分析 皮爾森相關(guān)分析結(jié)果表明，AK 含量與AP 含量、OM 含量、AN 含量的r值分別為0.363、0.173、0.194，這3 種土壤養(yǎng)分含量與AK 含量呈極顯著正相關(guān)，均為區(qū)域協(xié)同變量，在受外界環(huán)境影響方面具有相似性。AK 含量與E、地形起伏度（Topographic relief，TR） 、 地 表 粗 糙 度（Surface roughness，SR）這3 個(gè)地形要素的r值分別為-0.371、-0.221、-0.278，均呈極顯著負(fù)相關(guān)，證明地形為AK 含量的重要影響因子。另外，在研究區(qū)內(nèi)AK 含量隨海拔升高和地形起伏變大而降低。研究區(qū)位于三江平原，屬于溫帶大陸性氣候，降水多集中于夏季，對土壤侵蝕效應(yīng)明顯，受雨水及田塊微地形影響，AK 流失明顯。AK 含量與pH 值、坡度、坡向等無顯著相關(guān)關(guān)系，因此，在后續(xù)插值中剔除這些因素。

結(jié)合影像與實(shí)地調(diào)查發(fā)現(xiàn)，研究區(qū)部分耕層土壤長期深翻后白漿層被打破，在遙感影像上出現(xiàn)真彩色狀態(tài)下發(fā)白發(fā)亮的現(xiàn)象。通過NDVI、NDWI、RVI 等遙感指數(shù)與裸土?xí)r期影像特征對比分析發(fā)現(xiàn)，研究區(qū)內(nèi)裸土期NDWI、NDVI 與表層白漿化土壤在影像上的特征具有高度一致性。另將遙感指數(shù)與AK 含量進(jìn)行相關(guān)性分析后發(fā)現(xiàn)，裸土期NDWI、NDVI 與AK 含量的r值分別為0.164、-0.194，存在顯著相關(guān)關(guān)系，擇優(yōu)選取NDVI進(jìn)行COK插值。

2.2.2 典型白漿土區(qū)耕層土壤AK 輔助變量回歸分析 通過SPSS 26.0 進(jìn)行線性回歸分析，選取了AP、E、TR、NDVI等多種輔助變量進(jìn)行回歸分析，部分結(jié)果如表2 所示。結(jié)果顯示，上述輔助變量和研究區(qū)AK 含量存在顯著或極顯著相關(guān)關(guān)系。其中，AN 和AP的標(biāo)準(zhǔn)化系數(shù)分別為0.194、0.363，E、TR和SR標(biāo)準(zhǔn)化系數(shù)分別為-0.371、-0.221、-0.278，均與AK 含量呈極顯著相關(guān)；NDWI、NDVI 的標(biāo)準(zhǔn)化系數(shù)分別為0.164、-0.194，均與AK 含量呈顯著相關(guān)。該結(jié)果與皮爾森相關(guān)系數(shù)分析結(jié)果一致。

表2 典型白漿土區(qū)耕層土壤AK不同輔助變量回歸分析Tab.2 Regression analysis of different auxiliary variables of topsoil AK in typical albic soil area

2.2.3 典型白漿土區(qū)耕層土壤AK 半變異函數(shù)分析 在GS 9.0+中對AK 進(jìn)行半變異函數(shù)擬合分析發(fā)現(xiàn)，最優(yōu)擬合模型為球狀模型。半變異函數(shù)塊金值為0.04，變程為8.26 km。按照半方差分析法和空間相關(guān)性程度的分級標(biāo)準(zhǔn)，研究區(qū)較大范圍內(nèi)AK存在強(qiáng)烈空間變異。R2為0.886，表明模型具有良好的擬合優(yōu)度。

2.3 典型白漿土區(qū)耕層土壤AK模擬空間分析

為分析研究區(qū)內(nèi)AK的空間異質(zhì)性，將2種土壤養(yǎng)分因子（AP、AN）和2 種地形因子（E、TR）以及遙感指數(shù)（NDVI）這5 種因子作為輔助變量分別與AK進(jìn)行COK 插值。為了比較OK 插值與COK 插值的擬合精度，將AK 進(jìn)行OK 和COK 插值，得到6 幅AK空間分布圖，以可視化圖表征AK 空間異質(zhì)性，如圖3 所示。以AP、AN 和TR 為輔助變量的COK1、COK2和COK4插值結(jié)果均為斑塊狀分布，以E 和NDVI 為輔助變量的COK3和COK5插值結(jié)果均為條帶狀分布，且多種插值所得結(jié)果具有總體趨勢一致性，AK含量表現(xiàn)為南高北低。

圖3 典型白漿土區(qū)耕層土壤AK空間分布Fig.3 Spatial distribution of AK in typical albic soil area

交叉驗(yàn)證結(jié)果見表3。由多種不同驗(yàn)證數(shù)據(jù)來看，增加了輔助變量的克里格（協(xié)同克里格，COK）插值法優(yōu)于單一目標(biāo)變量的克里格（普通克里格，OK）插值。COK1—COK5的RMSE分別較OK 下降了2.096、0.992、0.300、1.459、0.173 mg/kg；MRE分別較OK 降 低 了1.51、0.25、0.22、0.65、0.08 個(gè) 百 分 點(diǎn)，MAE分別較OK 降低了0.50、0.16、0.07、0.35、0.02 mg/kg。從結(jié)果來看，以土壤養(yǎng)分因子（AP、AN）、地形因子（E、TR）和遙感指數(shù)（NDVI）作為輔助變量的COK 插值法精度更高，即COK＞OK，且以AP 作為輔助變量的效果最好。當(dāng)加入2種及以上因子作為輔助變量時(shí)COK 精度會(huì)有所降低，如將AP+TR 作為輔助變量時(shí)，COK 交叉驗(yàn)證RMSE為38.188 mg/kg，該結(jié)果精度低于基于AP 的COK（COK1，RMSE為36.225）mg/kg。以上表明，COK 插值并非輔助變量越多越好，應(yīng)該充分考慮與目標(biāo)信息的相關(guān)程度，擇優(yōu)選取協(xié)同變量，這樣才能獲得更高的擬合精度。

表3 典型白漿土區(qū)耕層土壤AK不同輔助變量插值精度系數(shù)Tab.3 Interpolation precision coefficient of different auxiliary variables of topsoil AK in typical albic soil area

3 結(jié)論與討論

一般認(rèn)為，土壤速效養(yǎng)分與農(nóng)業(yè)生產(chǎn)活動(dòng)密切相關(guān)[20]。曙光農(nóng)場內(nèi)AK 含量的空間變異系數(shù)為27.52%，與半變異函數(shù)的分析結(jié)果相一致，均屬于中等變異程度，并且與AK 在空間上的分布存在較強(qiáng)相關(guān)性。因此，認(rèn)為該地的AK 同時(shí)受自然環(huán)境因素和人為因素影響。2 種插值方式描述的AK 含量空間分布規(guī)律基本一致，均呈南高北低趨勢。插值結(jié)果除COK1、COK2和COK4插值均為斑塊狀分布外，其余均為條帶狀分布。從研究區(qū)整個(gè)農(nóng)場范圍來看，AK 的空間分布存在一定規(guī)律性，考慮到該地地形及耕地類型，研究區(qū)總體北高南低，且南部均為水田，各種土壤養(yǎng)分隨水分流失后聚集于南部區(qū)域，提高了南部區(qū)域多種養(yǎng)分含量。

地統(tǒng)計(jì)學(xué)中的克里格插值利用原始數(shù)據(jù)的空間依賴性，能夠?qū)崿F(xiàn)未知空間的無偏最優(yōu)估計(jì)，進(jìn)而達(dá)到土壤養(yǎng)分空間變異規(guī)律的最優(yōu)擬合，此方法多應(yīng)用于土壤屬性的空間異質(zhì)性研究[21]。與OK 插值法相比，COK 插值法將待測值無偏最優(yōu)估計(jì)方法從單一土壤屬性提高到2 個(gè)及以上協(xié)同變量，COK插值無論從模型擬合程度還是交叉驗(yàn)證精度來說都優(yōu)于OK插值[22]。

本研究結(jié)果顯示，無論輔助變量增加土壤養(yǎng)分因子（AP、AN）的COK1、COK2，還是增加地形因子（E、TR）的COK3、COK4，亦 或 是 增 加 遙 感 指 數(shù)（NDVI）的COK5，其RMSE、MRE和MAE值都低于OK。綜合比較后認(rèn)為，擬合精度表現(xiàn)為COK＞OK，證明了在同樣條件下添加具有較強(qiáng)相關(guān)性協(xié)同變量的COK 插值比OK 插值精度更高。這與趙業(yè)婷等[23]的研究結(jié)果相一致，COK 對土壤屬性變異特征的擬合優(yōu)度更高。本研究結(jié)果表明，COK 在農(nóng)場尺度典型白漿土區(qū)土壤養(yǎng)分屬性的模擬和預(yù)測中效果更優(yōu)。

在眾多環(huán)境變量中選取相關(guān)性強(qiáng)且更能反映目標(biāo)環(huán)境信息的輔助因子可獲取更加精確的土壤養(yǎng)分空間分布信息。本研究比較了添加土壤養(yǎng)分因子、地形因子、遙感指數(shù)與單一變量的OK 插值交叉精度驗(yàn)證結(jié)果，最終模型擬合效果為COK1＞COK4＞COK2＞COK3＞COK5＞OK。本研究選取E 作為輔助因子時(shí)插值效果不佳，可能是研究區(qū)位于三江平原地區(qū)，屬于典型白漿土區(qū)，土壤表層旱澇頻發(fā)，受降雨影響土壤被侵蝕后帶走大量養(yǎng)分，導(dǎo)致土壤本身原始養(yǎng)分含量降低，最終導(dǎo)致預(yù)測結(jié)果具有不穩(wěn)定性。雖然在本研究中坡度與AK 相關(guān)性較低，但以坡度為輔助變量的插值精度仍然較高。原因可能是數(shù)字高程模型分辨率不足以說明坡度與養(yǎng)分之間的關(guān)系。江葉楓等[24]的研究表明，輔助變量數(shù)量越多模擬精度越高。但本研究中以AP+TR 作為輔助變量的插值精度低于以AP 為輔助變量時(shí)的插值精度，證明了雖然輔助變量有助于模型精度的提高，但并非越多效果越好。與目標(biāo)變量相關(guān)性較小的輔助變量可能存在未知因素帶來的影響，甚至可能會(huì)對擬合精度產(chǎn)生負(fù)影響，這與馮家東等[19]的研究結(jié)果一致。另外，選取恰當(dāng)?shù)倪b感指數(shù)加入克里格插值也有助于提升插值精度。例如，加入NDVI 的COK5插值精度優(yōu)于單一的OK 插值精度。隨著遙感技術(shù)的發(fā)展，后續(xù)可進(jìn)一步分析其他遙感指數(shù)的適用性，以刻畫更加完善的環(huán)境變量，獲取更加詳盡的耕地土壤養(yǎng)分空間分布信息。

本研究主要利用地統(tǒng)計(jì)學(xué)方法并結(jié)合克里格空間插值技術(shù)，對典型白漿土區(qū)AK 的空間變異特征進(jìn)行分析。結(jié)果表明，曙光農(nóng)場耕層土壤AK 含量為92.00～364.00 mg/kg，平均為171.78 mg/kg，耕層土壤AK 含量整體處于較高水平?？臻g插值結(jié)果表明，研究區(qū)內(nèi)土壤鉀素含量南高北低，部分斑塊狀區(qū)域含量較高，為實(shí)現(xiàn)土壤改良，應(yīng)分區(qū)域調(diào)整鉀肥施用量。本研究中，模型擬合效果依次為COK1＞COK4＞COK2＞COK3＞COK5＞OK，以土壤養(yǎng)分因子AP為輔助變量的COK1較適合于典型白漿土區(qū)耕層AK的空間分布研究，而單一變量的OK 插值效果較差。綜上，選擇最優(yōu)輔助變量后的COK 插值法獲取典型白漿土區(qū)耕層土壤養(yǎng)分空間分布信息切實(shí)可行。