劉永坤,劉 偉,王昌陽(yáng),陳嘉明,馬 進(jìn)

(1. 江蘇師范大學(xué)地理測(cè)繪與城鄉(xiāng)規(guī)劃學(xué)院，江蘇 徐州 221116； 2. 中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221116； 3. 武漢大學(xué)遙感信息工程學(xué)院，湖北 武漢 430079)

土壤全鉀含量是植物生長(zhǎng)的主要營(yíng)養(yǎng)元素，同時(shí)也是土壤中影響農(nóng)作物產(chǎn)量的三要素之一[1]。研究土壤鉀含量的空間分布特征，不僅對(duì)均衡和擴(kuò)展土壤有效養(yǎng)分庫(kù)具有重要意義，還能為土地資源可持續(xù)利用和地區(qū)農(nóng)業(yè)的健康發(fā)展提供理論依據(jù)[2]。近幾年來(lái)，國(guó)外文獻(xiàn)[3—4]對(duì)土壤空間變異規(guī)律進(jìn)行了深層次的研究；國(guó)內(nèi)文獻(xiàn)[5]對(duì)土壤重金屬、土壤有機(jī)質(zhì)的空間變異進(jìn)行研究，取得了大量的研究成果。青海湖流域作為典型的復(fù)雜地貌類(lèi)型區(qū)，研究土壤全鉀含量對(duì)改善地區(qū)施肥效果、提高土壤利用效率具有一定的現(xiàn)實(shí)意義[6]。

土壤屬性的連續(xù)變化是土壤資源科學(xué)管理和使用的前提，土壤各類(lèi)屬性的空間插值是用來(lái)評(píng)估土壤性質(zhì)連續(xù)變化的主要方法[7]，同時(shí)也是“數(shù)字土壤”和“土壤計(jì)量學(xué)”的一個(gè)重要研究工具。空間插值精度取決于很多因素，一部分研究成果表明Kriging的模擬效果優(yōu)于反距離加權(quán)法(inverse distance weighting method,IDW)和Spline,而另一部分研究成果持相反觀(guān)點(diǎn)[8]。如文獻(xiàn)[9]通過(guò)比較3種插值方法評(píng)價(jià)土壤全氮含量的插值誤差，得出徑向基函數(shù)法插值模型優(yōu)于普通克里金插值法(ordinary kriging,OK)和IDW。修思玉等通過(guò)3種插值方法評(píng)價(jià)黑龍江省2010年32格氣象站點(diǎn)的年降雨量，采用平均誤差、平均絕對(duì)誤差、平均相對(duì)誤差等指標(biāo)進(jìn)行優(yōu)度檢驗(yàn)，得出Spline法最優(yōu)、IDW得到的降雨量的最大值與最小值與原始降雨量數(shù)據(jù)相差最大。文獻(xiàn)[10]對(duì)比了土壤Cd污染插值模型，得出不同的插值方法對(duì)極大值都存在不同程度的平滑。綜合以上研究發(fā)現(xiàn)，目前對(duì)空間屬性曲面建模的研究主要集中在單一全局插值模型。因此研究基于空間屬性的自適應(yīng)曲面建模具有重要的意義。

此外，一系列研究論證了利用地學(xué)環(huán)境變量為補(bǔ)充信息能夠有效提高插值精度。文獻(xiàn)[11]以高程、坡度、土地利用等輔助變量預(yù)測(cè)土壤精度。文獻(xiàn)[12]采用土壤電導(dǎo)率變量對(duì)不同尺度的土壤鹽分開(kāi)展了空間變異特征的隨機(jī)模擬。文獻(xiàn)[13]采用NDVI、高程和距河流距離為輔助變量研究黃河三角洲地區(qū)的鹽漬化，指出結(jié)合地學(xué)環(huán)境變量的插值模擬效果更好。但是，結(jié)合地學(xué)環(huán)境變量的插值模型也存在一些不足，如很多研究使用單一、同類(lèi)的地學(xué)環(huán)境變量解決不同地區(qū)的預(yù)測(cè)問(wèn)題，然而不同地區(qū)受到不同因素影響，其考慮的環(huán)境變量是不同的。各個(gè)插值方法的適應(yīng)性不同，不存在絕對(duì)最優(yōu)的空間插值模型，若能根據(jù)不同的插值模型自身的適應(yīng)性，分區(qū)進(jìn)行多模型集成，理論上可以提高模擬精度。

1 研究區(qū)概況與采樣方法

1.1 青海湖自然概況

研究區(qū)(36°38′N(xiāo)—37°29′N(xiāo),99°52′E—100°50′E)地處青海湖流域東南部。由于地殼運(yùn)動(dòng)形成了復(fù)雜多樣的地貌類(lèi)型，地貌以沖擊洪積扇平原為主，流域周?chē)笮『恿?0余條，屬于內(nèi)陸封閉水系，是許多野生動(dòng)物的棲息繁殖區(qū)域。青海湖流域?qū)俑咴敫珊档臏貛Т箨懶詺夂?，位于凍土地帶，干旱少雨，多風(fēng)，太陽(yáng)輻射強(qiáng)，溫差變化較大。研究區(qū)主要包括共和、剛察和海晏3個(gè)地區(qū)，除去青海湖，總面積約為2100 km2,海拔高度為3043～4516 m，有大量不同的土壤、植被、地貌類(lèi)型，屬于典型的復(fù)雜地貌類(lèi)型區(qū)。

1.2 樣品采集

研究通過(guò)空間分層組合的方式[14]進(jìn)行土壤采樣。在土壤采樣的同時(shí)記錄了采樣點(diǎn)的經(jīng)緯度、海拔、土地利用類(lèi)型、土壤類(lèi)型、植被類(lèi)型及地質(zhì)類(lèi)型等與土壤屬性密切相關(guān)的地學(xué)環(huán)境信息，每個(gè)樣點(diǎn)在0～15 cm、15～30 cm依次采樣兩次。筆者于2013年9月，在青海省環(huán)境檢測(cè)中心人員提供協(xié)助下，收集了青海湖典型地區(qū)土壤表層(0～30 cm)樣點(diǎn)數(shù)據(jù)110個(gè)。采樣結(jié)束后在實(shí)驗(yàn)室對(duì)樣品進(jìn)行風(fēng)干、研磨、過(guò)2 mm篩等流程后，取兩次土壤中全鉀元素的含量平均值作為記錄樣本值。

2 研究方法

2.1 地學(xué)環(huán)境要素篩選

已有研究表明，利用地學(xué)環(huán)境要素作為輔助變量能夠有效提高土壤全鉀含量插值精度和制圖效果。土壤全鉀空間變異性驅(qū)動(dòng)因子主要包括：土地利用類(lèi)型、土壤類(lèi)型、植被類(lèi)型、地質(zhì)類(lèi)型、坡度、施肥等[15-16]。根據(jù)已有的研究結(jié)論，結(jié)合本次主要針對(duì)自然景觀(guān)類(lèi)型區(qū)研究，排除施肥和土地管理措施2個(gè)驅(qū)動(dòng)因子，選取土地利用類(lèi)型、土壤類(lèi)型、植被類(lèi)型、地質(zhì)類(lèi)型4個(gè)地學(xué)環(huán)境要素作為輔助變量。其中土地利用類(lèi)型有灌叢、耕地和草原等6種；土壤類(lèi)型有栗鈣土、流動(dòng)和風(fēng)沙土等5種；植被類(lèi)型有水柏枝草場(chǎng)型、紫花針茅雜草草場(chǎng)型等30種；地質(zhì)類(lèi)型有沙丘、沙山和河谷平原等13種。使用SPSS計(jì)算輔助變量的土壤全鉀含量統(tǒng)計(jì)特征。為了進(jìn)一步篩選對(duì)土壤全鉀空間分布具有顯著影響的地學(xué)要素，利用SPSS軟件對(duì)土壤全鉀含量與上述4個(gè)地學(xué)要素進(jìn)行方差分析，選取具有顯著特征的地學(xué)要素作為輔助環(huán)境變量。表1的方差結(jié)果分析表明：在本文研究區(qū)域內(nèi)，植被類(lèi)型的景觀(guān)破碎度太高，很多區(qū)域內(nèi)只有1個(gè)采樣點(diǎn)，一些亞類(lèi)型的草地只有1～2個(gè)采樣點(diǎn)，模擬效果不理想，因此將植被類(lèi)型排除在外。

2.2 土壤屬性建模方法

2.2.1 插值方法的選取

本文中對(duì)土壤鉀含量進(jìn)行空間插值的4種方法均基于A(yíng)rcGIS 10.2平臺(tái)實(shí)現(xiàn)。傳統(tǒng)的插值方法有OK、IDW、回歸克里金(regression Kriging,RK)、貝葉斯克里金(Bayesian Kriging,BK)。OK在滿(mǎn)足二階平穩(wěn)假設(shè)和內(nèi)蘊(yùn)條件的基礎(chǔ)上，根據(jù)區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)(或協(xié)方差)的特點(diǎn)，確定待估點(diǎn)周?chē)阎c(diǎn)參數(shù)對(duì)待估點(diǎn)的加權(quán)值大小，再對(duì)待估點(diǎn)做出最優(yōu)估計(jì)的方法，其優(yōu)點(diǎn)是顧及各樣本點(diǎn)的空間相關(guān)性，在得到模擬結(jié)果的同時(shí)，得出估算精度的方差；IDW使用一組采樣點(diǎn)的線(xiàn)性權(quán)重組合來(lái)確定輸出的柵格值，權(quán)重與插值點(diǎn)和樣本點(diǎn)之間的距離成反比，是一種加權(quán)移動(dòng)平均的方法；RK綜合考慮了影響鉀空間變異的多種環(huán)境因子，但由于研究區(qū)環(huán)境要素分布比較破碎，導(dǎo)致沒(méi)有足夠多的采樣點(diǎn)來(lái)估計(jì)相對(duì)準(zhǔn)確的半方差函數(shù)。對(duì)全鉀含量進(jìn)行半變異分析，得到插值時(shí)的各參數(shù)值和擬合模型，其中塊金值較小，表明自身隨機(jī)性誤差引起的變異性不大；S/N+S比值接近于1，N/S比值小于30%，表明全鉀含量具有較強(qiáng)的空間相關(guān)性，擬合的最優(yōu)模型為指數(shù)模型[15]。

表1 不同地學(xué)要素類(lèi)型之間土壤全鉀的方差分析

注：*0.05顯著水平；**0.01顯著水平。

2.2.2 融合地學(xué)環(huán)境要素的自適應(yīng)曲面建模方法

2.2.2.1 融合地學(xué)環(huán)境要素的基插值模型

根據(jù)變異和空間相關(guān)理論，任何變量的空間變化都可以用下述兩個(gè)主要成分之和來(lái)表示：與局部變化有關(guān)的殘差成分、與趨勢(shì)有關(guān)的結(jié)構(gòu)成分。相應(yīng)的土壤全鉀含量的空間分布模式可以表示為

S(xi,l,k,yj,l,k)=direction(Geox,y)+residual(xi,l,k,yj,l,k)

(1)

式中,S(xi,l,k,yj,l,k)為第l種地學(xué)要素的第k種類(lèi)型的土壤全鉀采樣點(diǎn)模擬值，其中(xi,yj)為采樣點(diǎn)坐標(biāo)，i和j分別表示坐標(biāo)的行和列；direction(Geox,y)是描述(xi,yj)處第l種地學(xué)要素的第k種類(lèi)型S的趨勢(shì)值，其中Geo(x,y)是描述(xi,yj)處與土壤全鉀密切相關(guān)的地學(xué)環(huán)境信息；residual(xi,l,k,yj,l,k)是描述(xi,yj)處S的第l種地學(xué)要素的第k種類(lèi)型中第i行、第j列柵格點(diǎn)土壤全鉀的殘差值。基于均值模型得到趨勢(shì)函數(shù)，擬合描述結(jié)構(gòu)成分的趨勢(shì)面，實(shí)現(xiàn)趨勢(shì)分離；并根據(jù)土壤厚度采樣點(diǎn)的空間特征，比較選擇最佳的空間相關(guān)插值算法對(duì)殘差進(jìn)一步處理?；谝陨侠碚撝R(shí)，具體試驗(yàn)方案如圖1所示。

2.2.2.2 利用集成學(xué)習(xí)實(shí)現(xiàn)自適應(yīng)分區(qū)

根據(jù)上述建立的基插值模型Modeli對(duì)整個(gè)研究區(qū)進(jìn)行模擬生成一系列土壤全鉀插值曲面，使用采樣點(diǎn)土壤全鉀實(shí)測(cè)值減去預(yù)測(cè)值得到模擬誤差，以?huà)呙婢€(xiàn)算法(算法步驟如圖2所示)為集成學(xué)習(xí)的分類(lèi)學(xué)習(xí)器，對(duì)插值曲面進(jìn)行掃描，選擇分類(lèi)正確率較高且差異較大的掃描線(xiàn)進(jìn)行集成，對(duì)各個(gè)插值曲面進(jìn)行多種類(lèi)的自適應(yīng)分區(qū)，得到各個(gè)插值模型適用的空間范圍。

“+”和“-”分別表示滿(mǎn)足插值精度要求的樣點(diǎn)和不滿(mǎn)足插值精度要求的樣點(diǎn)。在這個(gè)過(guò)程中，使用水平和垂直的掃描線(xiàn)作為分類(lèi)器，對(duì)土壤全鉀插值曲面進(jìn)行分區(qū)。

(1) 根據(jù)對(duì)插值曲面分區(qū)的正確率，得到一個(gè)新的樣本分布(即樣本中每個(gè)樣本點(diǎn)的權(quán)重分布)D2、一個(gè)子分類(lèi)器h1。其中畫(huà)圈的樣本表示被分錯(cuò)的，比較大的“+”表示對(duì)該樣本作了加權(quán)。

(2) 根據(jù)分區(qū)的正確率，得到一個(gè)新的樣本分布D3、一個(gè)子分類(lèi)器h2。弱分類(lèi)器h2中有3個(gè)“-”符號(hào)分類(lèi)錯(cuò)誤，得到分區(qū)錯(cuò)誤率ε2=0.21，h2應(yīng)分配的權(quán)重α2=0.66。

(3) 得到一個(gè)子分類(lèi)器h3。弱分類(lèi)器h3中有2個(gè)“+”符號(hào)和一個(gè)“-”符號(hào)分類(lèi)錯(cuò)誤，得到分區(qū)錯(cuò)誤率ε3=0.14，h3應(yīng)分配的權(quán)重α3=0.91。

(4) 整合所有子分類(lèi)器，得到最終的Hfinal。經(jīng)過(guò)以上三步即可抽取所有滿(mǎn)足精度閾值的分區(qū)，每個(gè)區(qū)域?qū)儆谀膫€(gè)類(lèi)別，由該區(qū)域所在分類(lèi)器的權(quán)值綜合決定。經(jīng)過(guò)集成和NO Data區(qū)域的填充，得到土壤屬性自適應(yīng)曲面建模方法(adaptive surface modeling for soil properties,ASM-SP)模擬的空間分布圖。

2.3 插值結(jié)果精度檢驗(yàn)

獨(dú)立驗(yàn)證模型存在檢驗(yàn)精度不夠的問(wèn)題，使用交叉驗(yàn)證能夠更優(yōu)評(píng)價(jià)插值方法的質(zhì)量[17-19]。在交叉驗(yàn)證中，針對(duì)110個(gè)監(jiān)測(cè)點(diǎn)中的每一個(gè)點(diǎn)位數(shù)據(jù)，假設(shè)其數(shù)據(jù)未知，基于其余109個(gè)點(diǎn)位的數(shù)據(jù)，使用不同插值方法計(jì)算出模擬值，再結(jié)合觀(guān)測(cè)值求算誤差。

交叉驗(yàn)證的評(píng)定指標(biāo)有平均估計(jì)誤差百分比(PAEE)、相對(duì)均方差(RMSE)、均方根預(yù)測(cè)誤差(RMSPE)及殘差分析等。本文采用最大誤差(max error)、最小誤差(min error)、平均誤差(ME)、均方根誤差(RMSE)。

3 模擬結(jié)果與分析

3.1 土壤全鉀含量的統(tǒng)計(jì)特征

從研究區(qū)110個(gè)表層土壤鉀含量的統(tǒng)計(jì)結(jié)果分析(如圖3所示)可以看出，土壤鉀含量的值介于1.406 5%～2.346 4%之間，平均值為1.958 8%，K-S檢驗(yàn)結(jié)果表明鉀含量總體上服從正態(tài)分布。土壤鉀含量的變異系數(shù)為10.99%，表明研究區(qū)土壤鉀含量的空間變異性處于中等水平[15]。

3.2 ASM-SP

ASM-SP由3個(gè)步驟組成：①建立一系列基插值模型，即OK-Landuse，OK-Soil，OK-Geology；②基插值曲面的自適應(yīng)分割；③基插值優(yōu)化組合。自適應(yīng)分區(qū)方法詳見(jiàn)2.2節(jié)。ASM-SP構(gòu)建過(guò)程如下：

3.2.1 基插值模型建立

(1) 計(jì)算每種地學(xué)要素的土壤全鉀均值并得到趨勢(shì)曲面direction(Geox,y)，根據(jù)土壤全鉀實(shí)測(cè)值，結(jié)合輔助變量，得到土壤全鉀平均值(如圖3所示)。

(2) 用實(shí)測(cè)值減去土壤全鉀含量均值，計(jì)算得到土壤全鉀的模擬殘差。對(duì)模擬殘差進(jìn)行OK插值得到殘差曲面residual(xi,l,k,yj,l,k)(如圖4所示)。

(3) 利用式(1)計(jì)算第l種地學(xué)要素的第k種類(lèi)型的土壤全鉀采樣點(diǎn)模擬值S(xi,l,k,yj,l,k)，即將上述的趨勢(shì)曲面和殘差曲面相加。這是基插值曲面的構(gòu)造過(guò)程。

3.2.2 基于線(xiàn)性?huà)呙杷惴ǖ淖赃m應(yīng)曲面分割

基于第2.2.2節(jié)所述的構(gòu)造誤差曲面的方法，獲得不同插值模型的誤差面(如圖5所示)，并確定每個(gè)插值模型的適用范圍。

3.2.3 插值曲面的整合

在柵格單元優(yōu)化的基礎(chǔ)上，選取具有最小誤差的柵格單元的插值結(jié)果作為最佳柵格單元進(jìn)行集成。圖6顯示了與不同插值模型相對(duì)應(yīng)的最優(yōu)分區(qū)。

3.3 不同方法預(yù)測(cè)精度分析

為了評(píng)價(jià)ASM-SP模擬土壤全鉀含量空間分布的精度，本文比較了5種插值方法的精度，即RK、BK、IDW、OK+Landuse和ASM-SP。OK、BK結(jié)合輔助變量Landuse，與ASM-SP進(jìn)行對(duì)比；RK方法結(jié)合Landuse擬合趨勢(shì)項(xiàng)。逐點(diǎn)交叉驗(yàn)證的誤差指標(biāo)(見(jiàn)表2),對(duì)于最大誤差值來(lái)說(shuō)數(shù)值小精度高，對(duì)于最小誤差值來(lái)說(shuō)數(shù)值小精度高。表明ASM-SP預(yù)測(cè)誤差范圍最小，預(yù)測(cè)值的靈敏度和反映極值的情況優(yōu)于其他4種方法。針對(duì)ME和RMSE這幾個(gè)值而言，ASM-SP分別達(dá)到-0.002 1、0.281 1，相對(duì)其他結(jié)果較優(yōu)，且精度AC值最優(yōu)，因此ASM-SP顯示出更好的性能，其ME值比傳統(tǒng)插值方法(即OK和IDW)的ME值更接近于0。這意味著結(jié)合輔助變量的插值更加具有無(wú)偏性。并且ASM-SP的精度也充分表明其回歸曲線(xiàn)可以更好模擬預(yù)測(cè)值與真實(shí)值的關(guān)系。綜上所述ASM-SP的插值方法最優(yōu)。

總的來(lái)說(shuō)，ASM-SP插值方法優(yōu)于其他方法的主要原因有：①該方法結(jié)合了輔助變量，因此更準(zhǔn)確地刻畫(huà)土壤全鉀隨地學(xué)環(huán)境要素變化而突出的邊界；②使用線(xiàn)性?huà)呙杷惴▽⒆顑?yōu)區(qū)分割出來(lái)，并按此方法將各個(gè)地學(xué)要素得到最優(yōu)區(qū)域整合，在很大程度上提高了插值的精度。

表2 5種方法插值有效性交叉驗(yàn)證結(jié)果

3.4 不同預(yù)測(cè)方法效果對(duì)比

本文對(duì)比了5種土壤全鉀含量的插值效果(如圖7所示)。從圖7中可以看出，RK、BK和OK刻畫(huà)的變異程度較真實(shí)值區(qū)間小，插值顯示出不同程度的弱“牛眼”效應(yīng)。IDW插值可以體現(xiàn)出土壤全鉀分布的總體格局，但具有較強(qiáng)的“牛眼”效應(yīng)，且插值精度較低。

ASM-SP能夠有效地刻畫(huà)土壤全鉀空間變異格局，生成一個(gè)適中的插值范圍(1.30～2.32)，且能體現(xiàn)土壤全鉀空間分異較多的局部變化細(xì)節(jié)。BK、RK和OK插值得到的土壤全鉀值和IDW插值無(wú)法描述土壤屬性隨地學(xué)環(huán)境變量變化導(dǎo)致的突變邊界信息。ASM-SP方法對(duì)地形復(fù)雜地貌類(lèi)型區(qū)土壤屬性的空間插值具有較強(qiáng)的適應(yīng)性，能更準(zhǔn)確地描述研究區(qū)土壤屬性的空間變異格局。

4 結(jié) 語(yǔ)

利用地學(xué)環(huán)境變量對(duì)插值結(jié)果進(jìn)行修正，能夠解決土壤屬性空間分異突變的問(wèn)題，尤其能夠解決在丘陵溝壑等復(fù)雜地貌區(qū)和不同地學(xué)要素銜接區(qū)產(chǎn)生較大突變的問(wèn)題。土壤屬性的空間分異在較短水平距離內(nèi)產(chǎn)生較大的變化，導(dǎo)致單一全局插值模型精度有限。因此本文提出了一種融合地學(xué)環(huán)境要素，將不同插值模型最優(yōu)區(qū)進(jìn)行自適應(yīng)建模(ASM-SP)方法。在復(fù)雜地貌類(lèi)型區(qū)，ASM-SP能夠更準(zhǔn)確地刻畫(huà)土壤屬性的空間變異，有效降低模擬誤差。此外ASM-SP結(jié)合輔助變量，使其模擬結(jié)果更符合地學(xué)規(guī)律，便于物理解釋土壤屬性的空間變異特征。

對(duì)比未使用輔助變量的空間插值模型(IDW、OK)，以及利用地形因子輔助插值的RK法的結(jié)果顯示，ASM-SP刻畫(huà)土壤全鉀含量更加符合研究區(qū)土壤屬性空間變異規(guī)律，地學(xué)要素邊緣的細(xì)節(jié)信息表現(xiàn)得更加明顯。ME和RMSE等精度評(píng)價(jià)指標(biāo)也更小，分別達(dá)到-0.002 1和0.281 1，顯示出比其他插值模型更高的模擬精度，尤其能準(zhǔn)確刻畫(huà)復(fù)雜地貌類(lèi)型區(qū)土壤全鉀空間變異隨周?chē)貙W(xué)環(huán)境要素變化而突出的邊界。ASM-SP是實(shí)現(xiàn)復(fù)雜地貌類(lèi)型區(qū)土壤全鉀含量的高精度模擬的一種新方法，為以后的土壤屬性制圖研究提供了新思路和有益借鑒。