陳昌華，陳錫云，徐 英

（北京師范大學(xué) 地理學(xué)與遙感科學(xué)學(xué)院，北京100875）

土壤屬性的空間變異性研究對(duì)于實(shí)現(xiàn)土壤資源的科學(xué)管理與高效利用至關(guān)重要。土壤屬性空間分布圖是土壤屬性空間變異最直觀的表達(dá)形式，然而結(jié)合土壤剖面數(shù)據(jù)與土壤類(lèi)型圖制作土壤屬性分布圖的傳統(tǒng)方法帶有較強(qiáng)主觀性，難以定量反映土壤屬性的空間變異?？臻g統(tǒng)計(jì)是基于數(shù)據(jù)所代表的屬性具有空間自相關(guān)特性，根據(jù)已知空間數(shù)據(jù)來(lái)估計(jì)未知空間數(shù)據(jù)的數(shù)學(xué)方法［1］。20世紀(jì)70年代后期美國(guó)率先利用地理信息系統(tǒng)（GIS）的空間統(tǒng)計(jì)分析功能，開(kāi)展了土壤普查、土壤屬性數(shù)字制圖及土壤變異性方面的研究，推動(dòng)了利用地統(tǒng)計(jì)學(xué)方法研究土壤屬性空間變異的快速發(fā)展，使得反距離加權(quán)、普通克里金等方法在該領(lǐng)域的應(yīng)用趨于成熟。近年來(lái)，關(guān)于地統(tǒng)計(jì)插值方法比較及混合插值（hybrid interpolation）新模型［2］的發(fā)展和驗(yàn)證等方面的研究表明，以目標(biāo)變量的顯著相關(guān)性環(huán)境因子作為輔助信息（auxiliary information）的混合插值對(duì)提高不同區(qū)域的土壤屬性插值精度具有良好效果［3］。常用的輔助信息有地形因子［4］、氣候因子［5］及基于遙感衛(wèi)星影像的土地利用類(lèi)型［6］、土壤類(lèi)型等多分類(lèi)變量［7］。因輔助信息相對(duì)容易獲得，在一定程度上降低了野外工作量、提高了空間插值的預(yù)測(cè)能力，使得利用有限野外采樣點(diǎn)數(shù)據(jù)獲得更高精度的土壤屬性空間分布圖成為可能。

我國(guó)在20世紀(jì)80年代以來(lái)也有不少土壤學(xué)學(xué)者開(kāi)始這方面的研究，相關(guān)文獻(xiàn)報(bào)道主要集中在小尺度下確定性內(nèi)插和地統(tǒng)計(jì)插值方法比較［8－9］、土壤屬性插值效果評(píng)價(jià)及未采樣區(qū)土壤屬性預(yù)測(cè)等方面。對(duì)于大尺度低密度采樣的土壤普查數(shù)據(jù)，由于采樣不均勻且樣點(diǎn)分布較為稀疏，土壤屬性空間變異的隨機(jī)性因素增加，使類(lèi)似方法的應(yīng)用受限。但從理論上講，進(jìn)行空間插值的樣點(diǎn)數(shù)據(jù)滿足以下兩點(diǎn)［10］即可：一是空間樣點(diǎn)之間具有一定的相關(guān)性，二是空間變異具有連續(xù)性。因此，為充分發(fā)揮大尺度低密度采樣下土壤普查數(shù)據(jù)的潛在價(jià)值，利用全國(guó)第二次土壤普查東北地區(qū)的土壤全氮含量（STNC）數(shù)據(jù)，比較不同插值方法下土壤屬性的空間插值效果，同時(shí)利用土壤普查數(shù)據(jù)體系中顯著相關(guān)的環(huán)境因子數(shù)據(jù)作為輔助信息，探索提高土壤普查數(shù)據(jù)空間插值精度的有效途徑。為盡量避免土壤普查數(shù)據(jù)的采樣設(shè)計(jì)不合理及數(shù)據(jù)內(nèi)在變異性，比較了不同樣本容量下插值結(jié)果以充分認(rèn)識(shí)不同插值方法的空間趨勢(shì)預(yù)測(cè)穩(wěn)定性，希望能為今后土壤普查工作的土壤屬性制圖提供指導(dǎo)。

1 材料與方法

1.1 研究區(qū)域和數(shù)據(jù)來(lái)源

研究區(qū)域包括黑龍江、吉林、遼寧和內(nèi)蒙古東部的赤峰市、通遼市、興安盟和呼倫貝爾市（圖1），位于115°52′—135°09′E，38°72′—53°55′N(xiāo) 之間，總面積約為1.24×106km2，是世界三大黑土帶之一和我國(guó)廣義的黑土分布區(qū)。土壤數(shù)據(jù)來(lái)源于中國(guó)第二次土壤普查資料（1979—1992年），排除典型剖面中缺少經(jīng)緯度坐標(biāo)和STNC值的記錄，共417個(gè)有效樣點(diǎn)（圖1）。每個(gè)樣點(diǎn)還整理了詳細(xì)的地理位置、海拔高度、年均溫、年降水量等環(huán)境特征數(shù)據(jù)和土種類(lèi)型、所屬土類(lèi)、剖面厚度、pH、土壤陽(yáng)離子交換量等土壤理化數(shù)據(jù)，形成完整的輔助信息因子數(shù)據(jù)庫(kù)。

圖1 研究區(qū)域和有效樣點(diǎn)分布

1.2 數(shù)據(jù)處理

1.2.1 特異值處理與統(tǒng)計(jì)分析 由于特異值的存在會(huì)造成插值變量空間連續(xù)表面中斷，采用“樣本平均值3倍的標(biāo)準(zhǔn)差”來(lái)界定特異值，在該區(qū)間以外的數(shù)據(jù)均定為特異值，并分別用正常的最大值和最小值代替特異值［11］。利用SPSS 18.0軟件對(duì)特異值處理后的STNC數(shù)據(jù)進(jìn)行單樣本K—S（Kolmogorov—Smirnov）檢驗(yàn)和描述統(tǒng)計(jì)分析，以檢驗(yàn)數(shù)據(jù)分布情況。最后基于ArcGIS 10.0軟件平臺(tái)形成土壤氮素空間屬性數(shù)據(jù)庫(kù)，并對(duì)STNC數(shù)據(jù)作對(duì)數(shù)轉(zhuǎn)換計(jì)算，以滿足空間插值需要。

1.2.2 空間插值 采用確定性插值和地統(tǒng)計(jì)插值方法中較為常用且具有代表性的反距離加權(quán)法（inverse distance weighting，IDW）、徑向基函數(shù)（radial basis functions，RBF）、普 通 克 里 金 （ordinary kriging，OK）和回歸克里金（regression kriging，RK）4種方法對(duì)STNC進(jìn)行空間插值，前3種方法的相關(guān)報(bào)道較多［9－11］，其理論依據(jù)這里不再贅述。回歸克里金（RK）方法是一種可使用連續(xù)變量又可使用類(lèi)別變量作為輔助信息，同時(shí)結(jié)合Kriging插值方法來(lái)預(yù)測(cè)土壤屬性的混合插值技術(shù)?？捎梅匠淌剑?）—（2）表示［12］：

式中：f（x）——地理要素空間預(yù)測(cè)趨勢(shì)項(xiàng)，用一個(gè)隨輔助信息平穩(wěn)變化的線性函數(shù)表示；a0，a1——系數(shù)項(xiàng)，使預(yù)測(cè)模型適宜任何偏離率；m（x）——影響空間預(yù)測(cè)趨勢(shì)的輔助信息因子；Z（x）——土壤屬性預(yù)測(cè)值；ε（x）——f（x）趨勢(shì)預(yù)測(cè)土壤屬性值的殘差，模擬空間分布趨勢(shì)及隨機(jī)因素導(dǎo)致的不確定性。

本文采用趨勢(shì)項(xiàng)f（x）值與殘差項(xiàng)的Kriging預(yù)測(cè)值之和作為RK方法的預(yù)測(cè)結(jié)果，即STNC與環(huán)境變量之間的線性回歸預(yù)測(cè)值和回歸殘差值的OK插值之和作為RK方法的預(yù)測(cè)值［13］。其中，參與線性回歸分析的輔助信息因子共8個(gè)：經(jīng)度（longitude，Lo）、緯度（latitude，La）、海拔（H）、年均溫（mean annul temperature，MAT）、年均降水量（mean annul precipitation，MAP）、土層厚度（d）、pH 值和陽(yáng)離子交換容量（CEC）。

參考目前研究較為成熟的縣域尺度空間插值模型及本研究中不同樣本容量下樣點(diǎn)尺度的變化，在ArcGIS 10.0和 GS＋Version 7.0中，對(duì)4種空間插值方法在大尺度低密度采樣下的函數(shù)模型和參數(shù)設(shè)定詳見(jiàn)表1。OK和RK方法中涉及到Kriging插值均首先在GS＋軟件中進(jìn)行相應(yīng)變量的半方差分析，再基于以上所得最優(yōu)半方差模型參數(shù)，在ArcGIS軟件平臺(tái)下完成空間插值。相關(guān)研究表明：通過(guò)軟件得到的各向同性半方差模型實(shí)際上是對(duì)各向異性半方差函數(shù)模型均化的結(jié)果［11］，其差別不會(huì)對(duì)Kriging插值產(chǎn)生明顯影響，因此在OK和RK半方差分析中均采用各向同性半方差模型。IDW和RBF兩種方法由于擬合模型參數(shù)調(diào)試相對(duì)方便，直接在ArcGIS平臺(tái)上完成。

表1 不同插值方法模型基本參數(shù)設(shè)定情況

1.2.3 樣本容量設(shè)定 為充分認(rèn)識(shí)土壤調(diào)查數(shù)據(jù)STNC的空間變異性特征，以110［14］為最小個(gè)數(shù)設(shè)計(jì)了7個(gè)不同容量的樣點(diǎn)系列。具體操作為：利用Arc－GIS中的Create Subsets功能，隨機(jī)抽取417條樣點(diǎn)記錄中80%的記錄作為樣本A1，抽取A1中80%的記錄作為樣本A2，從A2中抽取80%的記錄作為樣本A3，…，依次按比例抽取形成6個(gè)不同容量的樣本數(shù)據(jù)，再加上全部的417條樣點(diǎn)記錄，共形成7個(gè)樣本系列。由于STNC原始數(shù)據(jù)中部分參與回歸分析的環(huán)境因子數(shù)據(jù)缺失，RK方法初始樣本容量為337，并單獨(dú)記為樣本B，之后從A3開(kāi)始計(jì)算，盡量保證4種空間插值方法使用統(tǒng)一樣點(diǎn)系列的STNC值進(jìn)行空間趨勢(shì)預(yù)測(cè)。

1.2.4 插值結(jié)果精度檢驗(yàn)與比較 選擇交叉驗(yàn)證方法中的平均誤差（ME，mean error）、均方根誤差（RMSE，root mean square error）和預(yù)測(cè)值與觀測(cè)值之間的相關(guān)系數(shù)R作為不同空間插值模型精度檢驗(yàn)指標(biāo)。為了消除不同樣本容量下STNC平均值和標(biāo)準(zhǔn)差的波動(dòng)性，將不同樣本容量下計(jì)算的ME和RMSE分別除以各自的平均值（mean）和標(biāo)準(zhǔn)差（SD）作標(biāo)準(zhǔn)化處理。標(biāo)準(zhǔn)化后的ME／mean和RMSE／SD值越小，R值越大則說(shuō)明空間插值精度越高。同時(shí)參照相關(guān)文獻(xiàn)［15］引入相對(duì)提高指數(shù)（RI，Relative Improvement）來(lái)定量化表示RK方法對(duì)RMSE精度的

式中：RMSE′——標(biāo)準(zhǔn)化 RMSE，即對(duì)應(yīng)的 RMSE／SD值；i——IDW，RBF，OK方法。

RI越大則說(shuō)明RK方法的精度提高程度越高，即輔助信息的添加越有效。提高度，其計(jì)算公式為：

2 結(jié)果與討論

2.1 STNC統(tǒng)計(jì)特征

417條有效樣點(diǎn)數(shù)據(jù)的統(tǒng)計(jì)分析（表2）表明，STNC原始取值范圍為0.08～21.48g／kg，剔除特異值后的取值范圍為0.08～5.31g／kg。

剔出特異值前后STNC變異系數(shù)均大于80%，表現(xiàn)出較強(qiáng)的變異性，這與研究區(qū)域自然條件復(fù)雜，土壤類(lèi)型和土地利用類(lèi)型多樣有關(guān)。但剔除特異值降低了數(shù)據(jù)的偏度和峰度，且使偏度和峰度更為接近。剔除特異值后的STNC進(jìn)行單樣本K—S檢驗(yàn)表明，對(duì)數(shù)轉(zhuǎn)換之前STNC屬于偏正態(tài)分布（p＜0.001）；經(jīng)過(guò)對(duì)數(shù)轉(zhuǎn)換后STNC的Z值由對(duì)數(shù)轉(zhuǎn)換之前的4.12減小為1.03，且 P 值變?yōu)?.244（＞0.01）。說(shuō)明經(jīng)對(duì)數(shù)轉(zhuǎn)換后STNC數(shù)據(jù)呈正態(tài)分布，空間趨勢(shì)預(yù)測(cè)中需要對(duì)有效樣點(diǎn)的STNC數(shù)據(jù)進(jìn)行對(duì)數(shù)轉(zhuǎn)換，以保證空間插值數(shù)據(jù)滿足正態(tài)分布。不同樣點(diǎn)系列之間，隨樣本容量減小，樣點(diǎn)最大間距減小，而最小間距的變化幅度不大（表2）。樣本系列A5以后最大采樣間距變化較大，但空間插值范圍仍保持A1時(shí)的90%以上，對(duì)STNC空間變異性在大尺度上的分辨不會(huì)存在影響。不同樣本容量下STNC基本描述統(tǒng)計(jì)量波動(dòng)不大，偏度系數(shù)在2.57～2.89之間；峰度系數(shù)在7.86～9.73之間，說(shuō)明各樣本系列具有相似的統(tǒng)計(jì)特征，不同容量的樣點(diǎn)系列設(shè)計(jì)比較合理。空間插值過(guò)程中不再針對(duì)每個(gè)樣本進(jìn)行正態(tài)分布檢驗(yàn)，空間預(yù)測(cè)完成后再將對(duì)數(shù)轉(zhuǎn)換后的數(shù)據(jù)計(jì)算結(jié)果進(jìn)行還原。

表2 STNC數(shù)據(jù)描述性統(tǒng)計(jì)特征及不同樣本系列數(shù)據(jù)特征比較

2.2 輔助信息因子的分析與篩查

根據(jù)相關(guān)性分析和回歸分析結(jié)果得到STNC回歸預(yù)測(cè)值及回歸殘差值，是采用RK方法分離趨勢(shì)項(xiàng)和殘差項(xiàng)的重要計(jì)算過(guò)程。由于進(jìn)行逐步回歸時(shí)因子變量間的共線性會(huì)影響結(jié)果的真實(shí)性，首先對(duì)8個(gè)環(huán)境變量進(jìn)行共線性判別（表3），以容忍度（tolerance，T）大于0.1和因子膨脹系數(shù)（variance inflation factor，VIF）小于10作為因子非共線性的判斷標(biāo)準(zhǔn)［16］。

結(jié)果表明，8個(gè)因子之間均滿足非共線性，可以全部作為自變量參與到STNC的逐步線性回歸分析中。通過(guò)相關(guān)性分析得到STNC與Lo，La，H，CEC呈顯著正相關(guān)，與MAT，d，pH呈顯著負(fù)相關(guān)，而與 MAP相關(guān)性不顯著。這可能是由于人類(lèi)耕作過(guò)程中灌溉措施在一定程度上降低了降水對(duì)土壤養(yǎng)分的影響。

表3 STNC與影響因子的相關(guān)性及變量共線性診斷分析

通過(guò)逐步線性回歸，設(shè)定因子變量進(jìn)入方程的概率水平為0.05，剔除的概率水平為0.10，得到STNC預(yù)測(cè)的最優(yōu)因子為CEC，MAT，d和pH，最優(yōu)線性回歸方程式如公式（4）。

考慮到本研究空間尺度較大且自然因素和土地利用變化較為復(fù)雜，同時(shí)對(duì)比其他研究尺度下RK方法回歸方程的決定系數(shù)［3，15，17］，該回歸擬合方 程的決定系數(shù)（R2＝0.328）處于可接受范圍。以此計(jì)算得到的STNC回歸預(yù)測(cè)殘差值呈近似正態(tài)分布（圖2），故用OK方法進(jìn)行回歸預(yù)測(cè)殘差值的空間插值時(shí)不再做數(shù)據(jù)轉(zhuǎn)換處理。

圖2 STNC回歸預(yù)測(cè)值的殘差分布情況

2.3 不同樣本容量下STNC空間變異性分析

利用最優(yōu)半方差函數(shù)模型定量化描述土壤屬性空間變異性，是利用地統(tǒng)計(jì)方法進(jìn)行土壤屬性空間插值的基礎(chǔ)。通過(guò)比較不同模型擬合的R2值和RSS值大小，STNC除樣本A3滿足球狀模型（Sh）外（表4），其他樣本容量下均滿足高斯模型（Gs），說(shuō)明樣本容量變化對(duì)擬合空間結(jié)構(gòu)的半方差函數(shù)模型影響不大。在實(shí)際參數(shù)設(shè)定過(guò)程中，有效滯后距和步長(zhǎng)的乘積為最大采樣間距的1／2，由于樣本容量的減少對(duì)樣點(diǎn)最大間距的影響，有效滯后距有所下降。塊金效應(yīng)（Nugget／Sill，N／S）是耕作措施、種植制度等隨機(jī)因素引起的空間異質(zhì)性占總變異的比例，當(dāng)樣本容量≥171時(shí)，STNC的 N／S 值為31.81%～49.93%，表現(xiàn)出中等水平的空間自相關(guān)性，且隨著樣本容量的減少N／S值基本呈現(xiàn)上升趨勢(shì)，說(shuō)明隨機(jī)性因素影響比重上升造成了空間自相關(guān)性的減弱。當(dāng)樣本容量＜171時(shí)，N／S值均小于25%，此時(shí)半方差模型對(duì)空間變異性的辨識(shí)能力下降使隨機(jī)性因素被掩蓋，空間變異結(jié)構(gòu)發(fā)生變化，得到的空間分析結(jié)果難以確信。STNC回歸預(yù)測(cè)殘差值的N／S取值范圍基本落在STNC塊金效應(yīng)波動(dòng)范圍以內(nèi)，說(shuō)明回歸預(yù)測(cè)殘差值保留了STNC固有空間結(jié)構(gòu)。STNC變程為1 643.72～3 635.57km，當(dāng)樣本容量≥171時(shí)，STNC樣點(diǎn)最大間距均落在空間自相關(guān)范圍內(nèi)，此時(shí)各樣本容量下的平均變程為2 040.53km。不同樣本容量下STNC回歸預(yù)測(cè)殘差值的空間自相關(guān)范圍也均大于2 000km。

表4 STNC半方差模型及擬合參數(shù)

目前針對(duì)大尺度低密度采樣的土壤屬性空間插值的報(bào)道很少，對(duì)比東北地區(qū)已有的相關(guān)研究，在小流域尺度（992hm2）上292個(gè)樣點(diǎn)的STNC空間插值的N／S塊金效應(yīng)為89.2%，變程為1 899m［18］?？h域尺度上，張素梅等［3］對(duì)自采的吉林省農(nóng)安縣471個(gè)樣點(diǎn)的STNC空間插值的N／S值為36.7%，變程為109.3km。Liu等［19］利用吉林省德惠縣的354個(gè)土壤肥力普查樣點(diǎn)數(shù)據(jù)所作的SOCC空間插值的N／S值為16.1%，變程達(dá)到632km。區(qū)域尺度上，Zhang等［20］利用995個(gè)自采樣點(diǎn)對(duì)黑龍江和吉林兩省黑土區(qū)的STNC空間插值的N／S值為5.97%，變程為1 037.5km。上述結(jié)果不難發(fā)現(xiàn)，不同空間尺度下不同樣本容量的土壤碳或氮的塊金效應(yīng)較為一致地表現(xiàn)出中強(qiáng)度空間相關(guān)性，且變程隨著各研究尺度的增大而增大。本研究區(qū)域?yàn)檎麄€(gè)東北地區(qū)，研究數(shù)據(jù)的樣點(diǎn)間距遠(yuǎn)大于上述同類(lèi)研究，得到的空間自相關(guān)范圍也較大。

2.4 不同方法的空間插值結(jié)果

比較4種插值方法的空間趨勢(shì)預(yù)測(cè)結(jié)果（圖3），整體來(lái)看各空間插值方法均能較直觀的反映出該區(qū)STNC空間分布的基本態(tài)勢(shì)，STNC均呈現(xiàn)從東北向西南逐漸遞減規(guī)律。高值區(qū)（＞1.10g／kg）呈倒“V”形條帶狀貫穿于研究區(qū)域的東北部地區(qū)，包括黑龍江省大部分地區(qū)和內(nèi)蒙古東北部的呼倫貝爾市及三江平原地區(qū)。所對(duì)應(yīng)的基本是暗棕壤、漂灰土、黑鈣土林區(qū)及沼澤濕地等理論上均有利于有機(jī)質(zhì)積累的區(qū)域。西北部中值區(qū)（0.63～1.10g／kg）包括松嫩平原在內(nèi)的吉林省和內(nèi)蒙古的興安盟、赤峰市以及通遼市的西北部地區(qū)。該區(qū)域主要為農(nóng)業(yè)耕作區(qū)和草甸草原，其中受人類(lèi)活動(dòng)影響的松嫩平原多發(fā)育一些草甸白漿土、白漿化暗棕壤等過(guò)渡土壤類(lèi)型，STNC水平適中。低值區(qū)（＜0.63g／kg）則呈斑塊狀分布在通遼市的東南部和遼寧省的西北部沿海地區(qū)，是沙漠、旱地和水田區(qū)，土壤類(lèi)型以風(fēng)沙土、潮土和濱海鹽土等有機(jī)質(zhì)積累較弱的土壤類(lèi)型為主。此插值結(jié)果分布趨勢(shì)與Yang等［21］使用全國(guó)第二次土壤普查數(shù)據(jù)得到的土壤有機(jī)碳密度的空間插值結(jié)果在該區(qū)的分布態(tài)勢(shì)保持很好的一致性。

圖3 STNC不同空間插值方法結(jié)果比較

對(duì)比不同樣本容量下插值結(jié)果（A1，A4，A73個(gè)樣本系列，圖3）可以發(fā)現(xiàn)，隨著樣本容量降低，4種方法插值結(jié)果的中、低值區(qū)范圍擴(kuò)大，高值區(qū)縮小且有了明顯的帶狀化表現(xiàn)，尤其是采樣點(diǎn)比較稀疏的東北部三江平原地區(qū)、內(nèi)蒙古東北和黑龍江西北部，過(guò)分的帶狀化削弱了空間變化特征。在空間分布的細(xì)節(jié)表達(dá)上，相同樣本容量下IDW和RBF插值圖中出現(xiàn)了特異值斑塊或“牛眼”現(xiàn)象，造成空間趨勢(shì)表達(dá)不連續(xù)，OK和RK方法由于Kriging插值模型對(duì)空間特異性具有一定的平滑作用，STNC空間變化規(guī)律明顯。同一種空間插值方法下，隨著樣本容量的減少，整體變化規(guī)律一致但趨勢(shì)變化的連續(xù)性變差。對(duì)比RK方法下A1和A4兩個(gè)樣點(diǎn)系列的插值含量水平和斑塊面積變化，在樣本容量較少（A4）情況下RK方法在空間趨勢(shì)預(yù)測(cè)上仍較為穩(wěn)定并在西北部和東北部的樣點(diǎn)稀疏區(qū)體現(xiàn)出更多的細(xì)節(jié)變化，說(shuō)明輔助信息的加入對(duì)這些區(qū)域的空間預(yù)測(cè)具有一定幫助。

2.5 不同空間插值方法的精度分析

插值結(jié)果的精度比較是評(píng)價(jià)不同插值方法有效性與適用性的重要手段，不同插值方法對(duì)土壤屬性插值精度在不同地區(qū)相關(guān)研究中的表現(xiàn)不一，部分研究表明Kriging的效果優(yōu)于IDW［22］和RBF中的Spline函數(shù)模型［23］，也有研究表明IDW和Kriging插值方法均有較好的表現(xiàn)［8］。僅針對(duì)不同Kriging插值方法的研究多體現(xiàn)為RK方法或其他添加輔助信息的Kriging插值方法和模型在提高插值精度方面表現(xiàn)較好［24］。比較本研究中4種空間插值方法的標(biāo)準(zhǔn)化ME，標(biāo)準(zhǔn)化RMSE，R三個(gè)精度指標(biāo)可見(jiàn)（圖4），其空間預(yù)測(cè)效果順序?yàn)椋篟K＞OK＞RBF＞IDW。不同空間插值方法下3個(gè)精度指標(biāo)值隨著樣本容量的變化規(guī)律不明顯，計(jì)算不同樣本容量下各精度指標(biāo)平均值，得到RK方法的ME值同OK，RBF和IDW方法相比平均降低0.18，0.17和0.16，RMSE值平均分別降低了0.15，0.18和0.19，R 平均分別提高了0.40，0.42和0.43。

從插值精度的變化幅度來(lái)看，隨著樣本容量的變化，RK和OK方法對(duì)不同樣本系列表達(dá)效果較為穩(wěn)定。由上述分析可知，當(dāng)樣本容量＜171時(shí)，STNC空間變異結(jié)構(gòu)發(fā)生變化，此時(shí)A6，A7兩個(gè)樣本容量下反映出的精度變化水平也不可信。

圖4 不同方法STNC空間插值精度隨樣本容量的變化

通過(guò)計(jì)算各樣本容量的平均RI值（圖4），RK方法相對(duì)于其他空間插值方法將預(yù)測(cè)精度分別提高了16.15%（OK），18.50%（RBF）和19.40%（IDW），表明利用環(huán)境因子作為輔助信息使該區(qū)空間插值精度得到了有效提高。相關(guān)研究中，Simbahan［4］等利用地表反射率、相對(duì)高程及土壤電導(dǎo)率3個(gè)環(huán)境因子作為輔助性的RK方法提高了CS元素制圖精度的16%。Zhang［7］等利用水流功率指數(shù)、濕潤(rùn)指數(shù)及坡度3個(gè)環(huán)境變量作為輔助信息的RK方法對(duì)土壤有機(jī)質(zhì)空間插值提高了13.93%的精度。Sumfleth［15］等利用相對(duì)高程和河網(wǎng)高程2個(gè)環(huán)境變量作為輔助信息的RK方法對(duì)STNC空間預(yù)測(cè)精度提高了13%。本研究中使用的輔助信息是土壤普查工作中必須測(cè)定的項(xiàng)目或者是相應(yīng)土壤點(diǎn)位上配套的氣象站點(diǎn)數(shù)據(jù)，同上述研究相比，在輔助信息易獲得性以及與土壤屬性的時(shí)空匹配性上具有明顯的優(yōu)勢(shì)，這對(duì)于今后提高土壤普查數(shù)據(jù)的空間插值精度價(jià)值巨大。

3 結(jié)論

（1）我國(guó)第二次土壤普查期間東北地區(qū)STNC在0.08～21.48g／kg之間，剔除特異值前后STNC數(shù)值水平具有較大變異性。大尺度低密度采樣下STNC表現(xiàn)出中等水平空間自相關(guān)性，但變程大于同區(qū)域的小尺度研究結(jié)果。目前針對(duì)大尺度（＞1∶100 000）低密度采樣的土壤屬性空間插值方法的研究較少，本研究得到的STNC空間變異性特征只能與同區(qū)域小尺度上的相關(guān)研究進(jìn)行類(lèi)比分析，得到的變程值對(duì)于STNC空間自相關(guān)范圍的指示意義難以進(jìn)行有效評(píng)價(jià)，但空間插值方法是利用土壤普查數(shù)據(jù)進(jìn)行土壤屬性制圖的有效途徑，因此對(duì)于大尺度低密度采樣的土壤屬性空間插值研究需要進(jìn)一步的加強(qiáng)。

（2）添加了輔助信息的回歸Kriging（RK）是對(duì)經(jīng)度（Lo）、緯度（La）、海拔（H）、年均溫（MAT）、年降水量（MAP）、土層厚度（d）、pH 值、土壤陽(yáng)離子交換量（CEC）8個(gè)環(huán)境變量與STNC值進(jìn)行線性回歸擬合，將得到的回歸擬合殘差值運(yùn)用普通Kriging（OK）進(jìn)行OK空間插值模擬，由OK模擬后的殘差值與回歸擬合的預(yù)測(cè)值之和得到RK方法的空間預(yù)測(cè)值。共線性判別表明，8個(gè)因子之間均滿足非共線性，可以全部作為自變量參與到STNC的逐步線性回歸分析中；通過(guò)逐步線性回歸得到本研究RK方法的輔助信息為CEC，MAT，d和pH值，其中，STNC與CEC呈顯著正相關(guān)，與 MAT，d，pH 值呈顯著負(fù)相關(guān)關(guān)系。

（3）樣本容量變化對(duì)不同插值方法的STNC空間趨勢(shì)預(yù)測(cè)影響不明顯。但在細(xì)節(jié)表達(dá)上，相同樣本容量下OK和RK方法的STNC空間插值結(jié)果空間變化特征更鮮明。隨著樣本容量的減少，同一種空間插值方法的空間表達(dá)更簡(jiǎn)略。RK方法對(duì)無(wú)樣點(diǎn)區(qū)的空間趨勢(shì)預(yù)測(cè)體現(xiàn)出更多的細(xì)節(jié)變化。此外，在以后的土壤普查工作中，應(yīng)該著重加強(qiáng)本研究中氮含量插值結(jié)果的高、中、低值區(qū)的過(guò)渡地帶以及內(nèi)蒙古東北部的風(fēng)沙地、黑龍江東部與東北部的多年凍土區(qū)和吉林東南沿海的鹽漬地等生態(tài)敏感地區(qū)的布點(diǎn)設(shè)計(jì)。

（4）不同樣本容量下4種方法空間插值精度總體表現(xiàn)為：RK＞OK＞RBF＞IDW，大尺度低密度采樣下地統(tǒng)計(jì)插值方法要優(yōu)于確定性內(nèi)插。RK方法相對(duì)于其他空間插值方法將預(yù)測(cè)精度提高程度（RI）分別為 16.15%（OK），18.50%（RBF）和 19.40%（IDW）。土壤普查過(guò)程中所收集的環(huán)境因子作為輔助信息即可使該區(qū)土壤氮含量的空間插值精度得到有效地提高，而且在低樣本容量下RK方法在空間趨勢(shì)預(yù)測(cè)上仍較為穩(wěn)定，預(yù)示著輔助信息的添加對(duì)大尺度低密度采樣下STNC的空間表達(dá)具有較好的效果。

［1］ Lam N S.Spatial Interpolation［M］∥International Encyclopedia of Human Geography.Elsevier Science Ltd.，2009.

［2］ Takata Y，F(xiàn)unakawa S，Akshalov K，et al.Spatial prediction of soil organic matter in northern Kazakhstan based on topographic and vegetation information［J］.Soil Science and Plant Nutrition，2007，53（3）：289-299.

［3］ 張素梅，王宗明，張柏，等.利用地形和遙感數(shù)據(jù)預(yù)測(cè)土壤養(yǎng)分空間分布［J］.農(nóng)業(yè)工程學(xué)報(bào)，2010，26（5）：188-194.

［4］ Simbahan G C，Dobermann A，Goovaerts P，et al.Fineresolution mapping of soil organic carbon based on multivariate secondary data［J］.Geoderma，2006，132（3／4）：471-489.

［5］ McKenzie N J，Ryan P J.Spatial prediction of soil properties using environmental correlation［J］.Geoderma，1999，89（1／2）：67-94.

［6］ Bui E N，Moran C J.A strategy to fill gaps in soil survey over large spatial extents an example from the Murray-Darling basin of Australia［J］.Geoderma，2003，111（1／2）：21-44.

［7］ Zhang Shiwen，Huang Yuanfang，Shen Chonggan，et al.Spatial prediction of soil organic matter using terrain indices and categorical variables as auxiliary information［J］.Geoderma，2012，171（S）：35-43.

［8］ Schloeder C A，Zimmerman N E，Jacobs M J.Comparison of methods for interpolating soil properties using limited data［J］.Soil Science Society of America Journal，2001，65（2）：470-479.

［9］ 石小華，楊聯(lián)安，張蕾.土壤速效鉀養(yǎng)分含量空間插值方法比較研究［J］.水土保持學(xué)報(bào)，2006，20（2）：68-72.

［10］ Webster R.Soil resources and their assessment［J］.Phil.Trans.R.London：B，1997，352（1356）：963-973.

［11］ 劉付程，史學(xué)正，于東升，等.基于地統(tǒng)計(jì)學(xué)和GIS的太湖典型地區(qū)土壤屬性制圖研究：以土壤全氮制圖為例［J］.土壤學(xué)報(bào)，2004，41（1）：20-27.

［12］ Stacey K F，Lark R M，Whitmore A P，et al.Using a process model and regression Kriging to improve predictions of nitrous oxide emissions from soil［J］.Geo-derma，2006，135：107-117.

［13］ Odeh I O A，McBratney A B，Chittleborough D J.Further results on prediction of soil properties from terrain attributes：Heterotopic Cokriging and regression-Kriging［J］.Geoderma，1995，67（3-4）：215-226.

［14］ Webster R，Oliver M A.Sample adequately to estimate variograms of soil properties［J］.Journal of Soil Science，1992，43（1）：177-192.

［15］ Sumfleth K，Duttmann R.Prediction of soil property distribution in paddy soil landscapes using terrain data and satellite information as indicators［J］.Ecological Indicators，2008，8（5）：485-501.

［16］ Menard S W.Applied logistic regression analysis［M］.Thousand Oaks：CA：SAGE Publication，Inc.，1995.

［17］ Ziadat F M.Analyzing digital terrain attributes to predict soil attributes for a relatively area［J］.Soil Science Society of America Journal，2005，69（5）：1590-1599.

［18］ Wei Jianbing，Xiao Duning，Zhang Xingyi，et al.Topography and land use effects on the spatial variation of soil organic carbon：A case study in a tipical small watershed of the black soil region in Northeast China［J］.Euras Soil Sci.，2008，41（1）：39-47.

［19］ Liu Dianwei，Wang Zongming，Zhang Bai，et al.Spatial distribution of soil organic carbon and analysis of related factors in croplands of the black soil region，Northeast China［J］.Agriculture，Ecosystems and Environment，2006，113（1／4）：73-81.

［20］ Zhang Xingyi，Sui Yueyu，Zhang Xudong，et al.Spatial variability of nutrient properties in black soil of Northeast China［J］.Pedosphere，2007，17（1）：19-29.

［21］ Yang Yuanhe，Mohammat A，F(xiàn)eng Jianming，et al.Storage，patterns and environmental controls of soil organic carbon in China［J］.Biogeochemistry，2007，84（2）：131-141.

［22］ Kravchenko A，Bullock D G.A comparative study of interpolation methods for mapping soil properties［J］.Agronomy Journal，1999，91（3）：393-400.

［23］ Panagopoulos T，Jesus J，Antunes M D C，et al.A-nalysis of spatial interpolation for optimising management of a salinized field cultivated with lettuce［J］.European Journal of Agronomy，2006，24（1）：1-10.

［24］ Hengl T，Heuvelink G，Alfred S.A generic framework for spatial prediction of soil variables based on regression-Kriging［J］.Geoderma，2004，120（1／2）：75-93.