王幼奇，張 興，趙云鵬，包維斌，白一茹*

(1 寧夏大學(xué)資源環(huán)境學(xué)院，銀川 750021；2 旱區(qū)特色資源與環(huán)境治理教育部國(guó)際合作聯(lián)合實(shí)驗(yàn)室，銀川 750021)

土壤陽(yáng)離子交換量 (CEC) 直接反映土壤保蓄、供應(yīng)和緩沖陽(yáng)離子養(yǎng)分的能力，同時(shí)也影響著土壤有機(jī)質(zhì)、酸堿度和土壤結(jié)構(gòu)等性質(zhì)[1-2]。因此CEC 通常被作為土壤質(zhì)量評(píng)價(jià)指標(biāo)之一，同時(shí)也是土壤施肥和改良等的重要依據(jù)。砂田是干旱半干旱區(qū)農(nóng)民在長(zhǎng)期生產(chǎn)實(shí)踐中探索形成的一種獨(dú)有的保護(hù)性耕作方式，具有保溫增滲、抑制土壤蒸發(fā)和防止侵蝕等作用[3]。近年來(lái)隨著寧夏硒砂瓜產(chǎn)業(yè)的快速發(fā)展，砂田在發(fā)揮其巨大經(jīng)濟(jì)和生態(tài)效益的同時(shí)，也面臨著由于長(zhǎng)期耕作導(dǎo)致土壤肥力下降、質(zhì)量退化等一系列嚴(yán)重問(wèn)題[4]。因此精準(zhǔn)預(yù)測(cè)砂田土壤CEC 空間分布規(guī)律對(duì)于防止砂田退化和提高土地生產(chǎn)力具有積極作用。

目前，基于地統(tǒng)計(jì)學(xué)的克里格插值方法在土壤要素空間預(yù)測(cè)中最常用[5]。普通克里格插值法(OK)對(duì)于樣點(diǎn)數(shù)量和樣點(diǎn)本身的數(shù)據(jù)質(zhì)量具有一定的依賴性，土壤自身具有的強(qiáng)變異性使得在復(fù)雜環(huán)境下使用OK法對(duì)土壤屬性進(jìn)行空間插值已不能滿足當(dāng)前要素空間預(yù)測(cè)對(duì)于精度的要求[6]。因此利用輔助信息協(xié)助變量進(jìn)行空間插值的方法如回歸克里格(RK)等得到了極大的應(yīng)用。RK 可以綜合多個(gè)變量對(duì)土壤屬性進(jìn)行空間插值，有效提高了預(yù)測(cè)精度，但土壤屬性的強(qiáng)空間變異性使得最小二乘法(OLS)這一全局模型無(wú)法很好地體現(xiàn)土壤屬性的局部特征[7]。地理加權(quán)回歸模型(GWR)作為一種有效處理回歸分析中空間非平穩(wěn)性和空間依賴性的局部模型，近年來(lái)被廣泛應(yīng)用到科學(xué)研究中[8-9]。如江振藍(lán)等[10]探討了GWR 模型在土壤重金屬高光譜預(yù)測(cè)中的適用性及局限性；王合玲等[11]應(yīng)用GWR模型對(duì)艾比湖土壤有機(jī)質(zhì)和土壤因子響應(yīng)關(guān)系的空間非平穩(wěn)性進(jìn)行了研究；袁玉蕓等[12]利用傳統(tǒng)回歸和GWR模型分析于田綠洲表層土壤鹽分及其影響因素的空間非平穩(wěn)性。以上研究均表明GWR模型對(duì)存在空間非平穩(wěn)性的數(shù)據(jù)具有更強(qiáng)的解釋能力和估計(jì)精度。地理加權(quán)回歸克里格法(GWRK)是將GWR 模型與RK 法相結(jié)合，綜合多個(gè)變量對(duì)土壤屬性進(jìn)行局部空間插值，提高了插值精度，能更好地反映出土壤屬性的局部變異情況[13]。目前，關(guān)于寧夏砂田土壤屬性空間插值的研究多使用普通克里格插值法[14-15]，對(duì)于將GWR 模型應(yīng)用在土壤要素空間插值上的研究較少。本文以砂田土壤為研究對(duì)象，通過(guò)分析土壤CEC 及其影響因素之間的相關(guān)性，驗(yàn)證和比較OK、RK 和GWRK 的空間制圖效果和插值精度，以期為土壤CEC 空間變異研究和土壤肥力管理提供科學(xué)依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于寧夏中衛(wèi)市香山鄉(xiāng)(圖1)，屬于寧夏中部干旱帶，地處 36°57′ ～ 37°07′N，104°56′ ～105°15′E，區(qū)內(nèi)氣候?yàn)闇貛Т箨懶詺夂颍嗄昶骄鶜鉁?.8 ℃，年均降水量247.4 mm，降雨多集中于7—9月。全區(qū)地勢(shì)南北高中間低，平均海拔約1 740 m，土壤類型主要為淡灰鈣土。天然植被主要以旱生灌木、半灌木、耐寒的蒿屬和禾本科草類為主。土壤基本理化性質(zhì)見(jiàn)表1。

圖1 研究區(qū)樣點(diǎn)分布Fig. 1 Distribution of soil sampling sites in study area

表1 研究區(qū)土壤基本理化性質(zhì)Table 1 Soil physicochemical properties in study area

1.2 樣點(diǎn)布設(shè)與數(shù)據(jù)采集

以香山地區(qū)行政區(qū)劃圖為底圖，在研究區(qū)域采用1.5 km × 1.5 km 網(wǎng)格布點(diǎn)方式進(jìn)行采樣。實(shí)際采樣過(guò)程中有些樣點(diǎn)落在村莊、道路等地，通常在附近200 m 內(nèi)進(jìn)行調(diào)整，并用GPS 記錄調(diào)整后的坐標(biāo)，共布設(shè)108 個(gè)樣點(diǎn)。使用土鉆采集根層土壤樣品(0 ～20 cm)。樣品帶回實(shí)驗(yàn)室后經(jīng)自然風(fēng)干，剔除植物殘根及石礫等，碾磨分別過(guò)10 目和60 目篩待測(cè)。CEC的測(cè)定采用乙酸鈉-火焰光度法[1]。土壤有機(jī)質(zhì)采用重鉻酸鉀氧化-外加熱法測(cè)定[2]。土壤機(jī)械組成采用激光粒度儀(MS3000，Malvern instruments)測(cè)定，粒徑分級(jí)標(biāo)準(zhǔn)采用美國(guó)制[16]。

1.3 研究方法

1.3.1 回歸克里格法 回歸克里格法(RK)考慮到土壤屬性空間變異驅(qū)動(dòng)因子的復(fù)雜性，將線性回歸與克里格插值相結(jié)合。當(dāng)目標(biāo)變量與輔助變量存在相關(guān)關(guān)系時(shí)，先通過(guò)對(duì)目標(biāo)變量和輔助變量的相關(guān)性分析和線性逐步回歸擬合，建立目標(biāo)變量與輔助變量的多元(或一元)回歸關(guān)系，得到趨勢(shì)項(xiàng)。然后對(duì)所得殘差項(xiàng)進(jìn)行半方差分析，并使用普通克里格法對(duì)殘差項(xiàng)進(jìn)行空間插值，最后使用柵格計(jì)算器將趨勢(shì)項(xiàng)與殘差項(xiàng)兩項(xiàng)相加，即為回歸克里格的插值結(jié)果。公式表達(dá)為：

式中： z(s) 為目標(biāo)變量在s處的預(yù)測(cè)值，m(s)為使用線性逐步回歸得到的趨勢(shì)項(xiàng)， ε(s)為使用普通克里格插值得到的殘差項(xiàng)。

1.3.2 地理加權(quán)回歸克里格法 地理加權(quán)回歸克里格法(GWRK)是通過(guò)對(duì)目標(biāo)變量和輔助變量進(jìn)行地理加權(quán)回歸擬合，得到局部回歸的殘差項(xiàng)，然后使用普通克里格插值法(OK)對(duì)所得殘差項(xiàng)進(jìn)行插值。GWRK 將RK 中的全局回歸(式2)轉(zhuǎn)換成地理加權(quán)回歸模型(GWR)的局部回歸(式3)，能更好地體現(xiàn)土壤屬性空間變異的局部變化。

式中：yi為樣點(diǎn)i的因變量；xik為第i個(gè)樣點(diǎn)上第k個(gè)變量的觀測(cè)值；(ui,vi)為樣點(diǎn)i的地理空間坐標(biāo)；β0為回歸的常數(shù)項(xiàng)；βk(ui,vi) 為第i個(gè)采樣點(diǎn)上的第k個(gè)回歸參數(shù)；εi為殘差項(xiàng)。如果βk(u,v)在空間保持不變，則模型(式3)就變?yōu)槿帜Ｐ?式2)。

1.3.3 模型精度檢驗(yàn)指標(biāo) 為了評(píng)價(jià)模型的預(yù)測(cè)精度，選取以下3 個(gè)指標(biāo)對(duì)模型進(jìn)行精度評(píng)價(jià)。平均誤差表示預(yù)測(cè)值與實(shí)測(cè)值偏差的算術(shù)平均值；均方根誤差表示預(yù)測(cè)值與實(shí)測(cè)值偏差的平方和觀測(cè)次數(shù)比值的平方根。均方根誤差對(duì)預(yù)測(cè)中特大或特小誤差反應(yīng)比較敏感；相對(duì)精度改進(jìn)值(RI)是表示衡量模型模擬效果是否優(yōu)于僅對(duì)實(shí)測(cè)值取平均值的指標(biāo)。

式中：ROK、RRK和RGWRK分別代表OK、RK 和GWRK擬合值與實(shí)測(cè)值的相關(guān)系數(shù)，RI 值為正值即表示RK、GWRK 較OK 的預(yù)測(cè)精度高，值越大說(shuō)明預(yù)測(cè)精度提高地越多；RI 值為負(fù)值則表示RK、GWRK預(yù)測(cè)精度低于OK。

1.4 數(shù)據(jù)處理與分析 采用SPSS 18.0 和Excel 對(duì)數(shù)據(jù)進(jìn)行描述統(tǒng)計(jì)分析、相關(guān)性分析和逐步回歸分析，使用GS+9.0 和ArcGIS 10.2 軟件對(duì)數(shù)據(jù)進(jìn)行地統(tǒng)計(jì)分析和空間插值。為排除異常值對(duì)半方差函數(shù)穩(wěn)健性的影響，依據(jù)3σ 準(zhǔn)則對(duì)原始取樣數(shù)據(jù)進(jìn)行了異常值識(shí)別，異常值用正常的最大或最小值代替，后續(xù)分析采用處理過(guò)的原始數(shù)據(jù)進(jìn)行計(jì)算。為驗(yàn)證模型的預(yù)測(cè)精度，隨機(jī)選擇86 個(gè)樣點(diǎn)作為建模樣點(diǎn)進(jìn)行空間插值，22 個(gè)樣點(diǎn)作為驗(yàn)證樣點(diǎn)用于分析插值精度。

2 結(jié)果與分析

2.1 土壤CEC 描述統(tǒng)計(jì)分析

描述統(tǒng)計(jì)表明(表2)，建模點(diǎn)和驗(yàn)證點(diǎn)CEC 均值分別為9.82、10.47 cmol/kg，含量變化范圍分別為5.75 ～ 14.87 cmol/kg 和6.15 ～ 15.83 cmol/kg。研究區(qū)土壤CEC 均值為10.145 cmol/kg。根據(jù)土壤保肥能力分級(jí)方法[17]，供肥保肥能力弱的樣本數(shù)占總樣本數(shù)的60.19%，供肥保肥能力中等的占39.81%，說(shuō)明研究區(qū)土壤供肥保肥能力多處于低水平。偏度和峰度數(shù)值均較小(＜1)。從變異程度看，建模點(diǎn)和驗(yàn)證點(diǎn)CEC變異系數(shù)分別為17.57% 和21.02%，屬中等變異。經(jīng)K-S 非參數(shù)檢驗(yàn)，建模點(diǎn)和驗(yàn)證點(diǎn)的P值均大于0.05，數(shù)據(jù)服從正態(tài)分布。

2.2 土壤CEC 與土壤理化性質(zhì)的相關(guān)性分析

影響土壤CEC 含量的因素有很多，有機(jī)質(zhì)及礦質(zhì)膠體的數(shù)量與性質(zhì)是決定土壤CEC 的重要因子[2]。相關(guān)學(xué)者亦研究發(fā)現(xiàn)CEC 與有機(jī)質(zhì)、土壤質(zhì)地存在顯著相關(guān)性，明顯影響CEC 在區(qū)域上的分布特征[18-19]。因此為提高CEC 預(yù)測(cè)精度、更好展現(xiàn)其空間分異特征，本文從CEC 與有機(jī)質(zhì)、礦質(zhì)膠體的耦合關(guān)系為出發(fā)點(diǎn)，選取了研究區(qū)土壤有機(jī)質(zhì)、質(zhì)地與CEC 進(jìn)行相關(guān)性分析[20]。表3 可以看出，CEC 與有機(jī)質(zhì)、砂粒、粉粒和黏粒呈顯著相關(guān)關(guān)系。其中CEC 與有機(jī)質(zhì)、黏粒和粉粒呈顯著正相關(guān)，與砂粒呈顯著負(fù)相關(guān)。

表2 土壤CEC 的描述統(tǒng)計(jì)特征Table 2 Descriptive statistics of soil CEC

2.3 土壤CEC 的回歸模型與半方差模型

2.3.1 逐步回歸分析 在相關(guān)分析基礎(chǔ)上，選擇有機(jī)質(zhì)、黏粒、粉粒和砂粒作為輔助變量，采用逐步回歸篩選出對(duì)各個(gè)回歸模型具有顯著貢獻(xiàn)的回歸參數(shù)。逐步回歸結(jié)果表明(表4)，有機(jī)質(zhì)和砂粒是對(duì)CEC 插值的最佳自變量。方差膨脹因子(VIF)為1.12，小于7.5，說(shuō)明自變量間不存在共線性。調(diào)整后的決定系數(shù)為0.567，說(shuō)明OLS 模型對(duì)CEC 的方差解釋度為56.7%。擬合模型如下：y= 7.478+0.964x1-0.052x2，式中，y為CEC；x1為有機(jī)質(zhì)含量；x2為砂粒含量。為了便于對(duì)比，應(yīng)用GWR 模型時(shí)，同樣選用有機(jī)質(zhì)和砂粒進(jìn)行建模。

表3 土壤CEC 與土壤理化性質(zhì)的相關(guān)性Table 3 Correlation between soil CEC and other soil physicochemical properties

表4 土壤CEC 與土壤理化性質(zhì)的逐步回歸過(guò)程Table 4 Stepwise linear regression process of soil CEC with other physicochemical properties

2.3.2 地統(tǒng)計(jì)分析 對(duì)CECOLS 殘差和GWR 殘差進(jìn)行描述性統(tǒng)計(jì)分析。OLS 殘差變化范圍為 -1.49 ～4.02，均值為0.004 7；GWR 殘差變化范圍為 -1.57 ～4.44，均值為 -0.0345。經(jīng)K-S 非參數(shù)檢驗(yàn)，兩種回歸方法所得殘差均服從正態(tài)分布。在GS+9.0 中對(duì)實(shí)測(cè)值、OLS 殘差和GWR 殘差進(jìn)行半方差函數(shù)模型擬合(表5)?？梢钥闯?，三者均可用球狀模型擬合。塊金值(C0)代表隨機(jī)因素帶來(lái)的變異，通常由區(qū)域變量的變異性和測(cè)定誤差所造成[21]。三者的塊金值(C0)均小于0.3，表明CEC 由測(cè)定誤差或土壤性質(zhì)帶來(lái)的隨機(jī)變異較小，空間變異主要受結(jié)構(gòu)性因子影響。塊金系數(shù)(C0/(C0+C))表示由隨機(jī)因素引起的空間變異占系統(tǒng)總變異的比例。一般來(lái)說(shuō)，比值＜25% 時(shí)空間變異性表現(xiàn)為強(qiáng)空間相關(guān)性，介于25% ～ 75% 時(shí)表現(xiàn)為中等強(qiáng)度空間相關(guān)性，比值＞75% 時(shí)表現(xiàn)為弱空間相關(guān)性[17]。實(shí)測(cè)值、OLS 殘差和GWR 殘差塊金系數(shù)分別為8.50%、6.36% 和7.02 %，比值均在25% 以下，說(shuō)明CEC 及其殘差具有強(qiáng)烈空間自相關(guān)。

表5 土壤CEC、OLS 殘差和GWR 殘差的半方差模型參數(shù)Table 5 Semivariance parameters of soil CEC and regression residuals of OLS and GWR

2.4 土壤CEC 的空間分布特征

利用OK、RK 和GWRK 3 種方法進(jìn)行空間插值得到的CEC 空間分布圖。根據(jù)圖2 可以看出，3 種方法得到的結(jié)果中CEC 的分布格局大致相似，即研究區(qū)西部的CEC 要大于東部，高值區(qū)分布在西北、東北和南部，CEC 存在從這些地區(qū)向中東部逐漸減小的趨勢(shì)。CEC 預(yù)測(cè)值范圍分別為5.75 ～ 14.88、5.77 ～ 13.00、5.52 ～ 12.90 cmol/kg，各方法基本一致。從制圖效果看，OK 法插值圖較為粗略，RK 法的插值結(jié)果體現(xiàn)了CEC 隨著有機(jī)質(zhì)和砂粒的變化情況，但是斑塊的邊界較為零碎，這與現(xiàn)實(shí)中CEC 的空間漸變特征不相符。GWRK 法不僅保留了回歸分析中CEC 隨土壤理化性質(zhì)變化的豐富細(xì)節(jié)，還反映出CEC 的空間漸變特征，斑塊邊界更為光滑，感官上更符合CEC 的實(shí)際分布狀況。為了更好地反映3 種插值模型對(duì)區(qū)域土壤CEC 的預(yù)測(cè)精度，利用平均誤差、均方根誤差和RI 值對(duì)插值結(jié)果進(jìn)行精度評(píng)價(jià)。

2.5 不同插值方法的空間預(yù)測(cè)精度評(píng)價(jià)

表6 對(duì)OK、RK 和GWRK 方法進(jìn)行了精度評(píng)價(jià)。結(jié)果表明3 種方法的平均誤差基本接近于0，表明了所構(gòu)建模型的預(yù)測(cè)具有無(wú)偏性；均方根誤差OK＞RK＞GWRK，范圍為1.084 ～ 1.853，趨近于1，說(shuō)明所構(gòu)建模型的誤差有效。從相對(duì)精度改進(jìn)值來(lái)看，RK 和GWRK 法相較于OK 法擬合精度分別提高了40.49%和41.50%，表明3 種插值方法中GWRK法插值精度最高，為最優(yōu)插值模型。

圖2 三種插值方法下土壤CEC 的空間分布Fig. 2 Spatial distribution of soil CEC by OK, RK and GWRK interpolation

表6 OK、RK 和GWRK 方法的精度評(píng)價(jià)指標(biāo)對(duì)比Table 6 Precision evaluation of three interpolation methods

3 討論

研究區(qū)土壤CEC 均值為10.145 cmol/kg，供肥保肥弱的樣本數(shù)占總樣本數(shù)的60.19%，研究區(qū)土壤保肥供肥能力較弱。前人的相關(guān)研究也發(fā)現(xiàn)了類似情況[4,15]。其產(chǎn)生原因可能與當(dāng)?shù)赝寥滥纲|(zhì)和耕作制度有關(guān)。研究區(qū)處于我國(guó)西北內(nèi)陸腹地，土壤類型為低肥力的淡灰鈣土，自然環(huán)境條件限制了當(dāng)?shù)赝寥赖姆柿λ?。再加上研究區(qū)有覆砂的傳統(tǒng)耕作制度，農(nóng)民施肥難度大，外源性的有機(jī)肥施入量不夠，導(dǎo)致當(dāng)?shù)赝寥婪柿ζ?。砂田土壤CEC 在空間上展現(xiàn)出高度變異性，在研究區(qū)地形、氣候和母質(zhì)相對(duì)均一的情況下，土壤理化性質(zhì)可能是導(dǎo)致土壤CEC 在空間上呈現(xiàn)出較強(qiáng)變異性的重要原因。其中有機(jī)質(zhì)中膠體成分在土壤固相中所占比例較大，土壤膠體是CEC 的載體，隨有機(jī)質(zhì)含量增加，CEC 也隨之增加。同時(shí)在土壤顆粒組成中，黏粒膠體集中了80% 以上的負(fù)電荷，土壤黏粒含量愈高，CEC 就越高[18-20]。相關(guān)學(xué)者亦研究發(fā)現(xiàn)CEC 與有機(jī)質(zhì)、土壤質(zhì)地存在顯著相關(guān)性，明顯影響CEC 在區(qū)域上的分布特征。砂田CEC 與土壤有機(jī)質(zhì)、黏粒和粉粒呈正相關(guān)，與砂粒呈負(fù)相關(guān)，即在有機(jī)質(zhì)和黏粒含量高的地方，CEC相對(duì)充足。針對(duì)這一情況，建議在后續(xù)農(nóng)業(yè)生產(chǎn)過(guò)程中增施有機(jī)肥，推廣輪作制度以培肥土壤，從而提高作物產(chǎn)量。

空間插值技術(shù)是利用采樣點(diǎn)數(shù)據(jù)對(duì)研究區(qū)內(nèi)未知點(diǎn)進(jìn)行預(yù)測(cè)的方法。常用的空間插值方法有普通克里格法、回歸克里格法等[22-23]。張慧智等[24]對(duì)比了普通克里格、泛克里格和回歸克里格在土壤溫度空間預(yù)測(cè)中的效果，發(fā)現(xiàn)回歸克里格法預(yù)測(cè)精度最高，能夠更好地表達(dá)復(fù)雜地形地區(qū)的局部變異。Hengl 等[25]采用回歸克里格對(duì)土壤有機(jī)質(zhì)含量進(jìn)行空間預(yù)測(cè)，發(fā)現(xiàn)其預(yù)測(cè)精度較普通克里格表現(xiàn)更好，均方根誤差數(shù)值更低。但是，土壤理化性質(zhì)與土壤CEC 之間的相關(guān)關(guān)系因地而異，局部區(qū)域的CEC 與土壤理化性質(zhì)關(guān)系可能與做全局分析時(shí)得出的結(jié)論相反，這在回歸克里格中沒(méi)有考慮。GWR 模型能將數(shù)據(jù)的空間信息納入分析過(guò)程，通過(guò)計(jì)算回歸模型的局部參數(shù)來(lái)解決地理數(shù)據(jù)中存在的空間自相關(guān)性及空間非平穩(wěn)性問(wèn)題，從而提高模型的擬合優(yōu)度及模擬效果[26]。本文對(duì)3 種插值模型的精度分析結(jié)果也證明了GWRK 法在空間制圖方面表現(xiàn)出更豐富的局部細(xì)節(jié)和更光滑的空間漸變特征。因此，在后續(xù)土壤屬性的空間制圖研究中，應(yīng)注重GWRK 法在土壤屬性空間插值中的應(yīng)用。

4 結(jié)論

1) 相關(guān)分析表明，土壤CEC 與有機(jī)質(zhì)、黏粒和粉粒含量呈顯著正相關(guān)，與砂粒含量呈顯著負(fù)相關(guān)。

2) 地統(tǒng)計(jì)分析表明，土壤CEC 及其插值殘差均可用球狀模型擬合；實(shí)測(cè)值、OLS 殘差和GWR 殘差塊金值均小于0.3，說(shuō)明其空間變異主要受結(jié)構(gòu)性因子影響；塊金系數(shù)分別為8.50%、6.36% 和7.02%，表現(xiàn)出強(qiáng)烈空間自相關(guān)。

3)從插值效果和精度分析來(lái)看，GWRK 是3 種模型中插值效果最為理想的；RK 和GWRK 法相較于OK 法擬合精度分別提高了40.49% 和41.50%，3種插值方法中GWRK 法插值精度最高，為最優(yōu)插值模型。