基于地形分區(qū)的黃土丘陵區(qū)土壤有機(jī)質(zhì)空間預(yù)測(cè)

張 婧，畢如田，丁皓希，文偉杰，荊穎薔

（山西農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院，山西太谷 030801）

土壤有機(jī)質(zhì)作為表征耕地肥力大小的重要指標(biāo)，其含量在農(nóng)作物的整個(gè)生長(zhǎng)過(guò)程中發(fā)揮著非常重要的作用[1]。土壤有機(jī)質(zhì)具有較強(qiáng)的時(shí)空變異性，其空間規(guī)律一直是許多學(xué)者研究的熱點(diǎn)問(wèn)題。大量研究表明，土壤有機(jī)質(zhì)的空間變異會(huì)受各種因素的影響，如地形、土壤類型、耕作模式、土地利用方式等[2-3]。但是關(guān)于多種因素協(xié)同對(duì)有機(jī)質(zhì)的空間影響目前還少有研究。

Kriging 空間插值法是一種土壤屬性空間預(yù)測(cè)最普遍的方法[4-5]，但是普通Kriging 插值在平坦、地形單一的區(qū)域預(yù)測(cè)結(jié)果較好，而在復(fù)雜地區(qū)插值效果一般。隨著科學(xué)技術(shù)的進(jìn)步，Kriging 插值引入其他參數(shù)作為輔助變量的方法已經(jīng)得到廣泛的應(yīng)用，而且預(yù)測(cè)精度較高。例如，徐占軍等[6]基于分區(qū)Kriging 插值法對(duì)煤炭開(kāi)采沉陷區(qū)土壤有機(jī)碳含量進(jìn)行空間預(yù)測(cè)，結(jié)果顯示，分區(qū)Kriging 預(yù)測(cè)結(jié)果的精度比普通Kriging 預(yù)測(cè)的更高；MISHRA 等[7]以印第安納州為研究區(qū)，利用普通Kriging 模型和剖面深度分布函數(shù)分別預(yù)測(cè)了表層土壤有機(jī)碳的含量，結(jié)果表明，剖面深度分布函數(shù)預(yù)測(cè)精度高于普通Kriging 模型。

地形是引起土壤有機(jī)質(zhì)空間變異的主要因子，而地形起伏度是表征地形特征的重要參考指標(biāo)之一[8]。地形起伏度是表示某一確定區(qū)域中最高海拔和最低海拔之差[9]，它是描述區(qū)域地形特征的一個(gè)宏觀性指標(biāo)，用以表征地面的起伏狀況和切割程度，被廣泛應(yīng)用于自然地質(zhì)災(zāi)害等定量評(píng)價(jià)方面[10-14]，目前，結(jié)合地形起伏度進(jìn)行有機(jī)質(zhì)空間規(guī)律的研究還相對(duì)較少。

本試驗(yàn)以黃土丘陵區(qū)為研究區(qū)，選取地形起伏度和地貌因子進(jìn)行地形分區(qū)，然后進(jìn)行Kriging 插值，并將其與直接插值進(jìn)行比較，對(duì)比2 種方法的預(yù)測(cè)精度，旨在提出一種更適合丘陵山區(qū)土壤有機(jī)質(zhì)的插值預(yù)測(cè)方法，為獲取更精準(zhǔn)的土壤有機(jī)質(zhì)空間分布提供指導(dǎo)性建議。

1 材料和方法

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

1.1.1 研究區(qū)概況 試驗(yàn)以晉中市太谷區(qū)為研究區(qū)，該區(qū)位于山西省中部，地處晉中盆地東北部（東經(jīng)112°28′～113°01′，北緯37°12′～37°03′），隸屬于黃土高原，地形復(fù)雜多樣，主要有平原、丘陵和山地，由東南向北傾斜，海拔高度為763～1 900 m；土地利用類型主要為耕地、園地、林地和草地（圖1-A）。該區(qū)總體格局為西北低東南高，平原與丘陵交錯(cuò)分布；全區(qū)耕地面積為28 885.46 hm2，空間分布比較明顯，主要分布在763～933 m 的海拔范圍內(nèi)（圖1-B）；成土母質(zhì)主要為黃土狀物質(zhì)、河流沖積物、馬蘭黃土等，土壤類型為褐土和潮土。

1.1.2 數(shù)據(jù)來(lái)源 本研究所需數(shù)據(jù)主要包括分辨率為30 m 的SRTMDEM數(shù)據(jù)，來(lái)源于中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云網(wǎng)站（www.gscloud.cn）；土壤有機(jī)質(zhì)數(shù)據(jù)，來(lái)源于山西省2010 年測(cè)土配方項(xiàng)目采樣；耕地?cái)?shù)據(jù)，來(lái)源于2015 年1∶10000 太谷區(qū)土地利用變更調(diào)查數(shù)據(jù)庫(kù)以及2017 年太谷區(qū)耕地質(zhì)量等級(jí)更新成果和山西地貌。

1.2 研究方法

1.2.1 地形分區(qū)

1.2.1.1 地形起伏度的提取 使用Arcgis 矩形窗口大小為n×n 像元（n＝2，3，…，23）提取地形起伏度：一統(tǒng)計(jì)n×n（n＝2，3，…，23）不同窗口內(nèi)像元的最大值和最小值；二計(jì)算不同窗口內(nèi)最大值與最小值之間的差值；三統(tǒng)計(jì)各窗口下的平均起伏度值。

1.2.1.2 最佳地形起伏單元的計(jì)算 均值變點(diǎn)分析法檢驗(yàn)最有效的是恰有一個(gè)變點(diǎn)[15]，已經(jīng)被廣泛應(yīng)用于提取地形起伏度最佳統(tǒng)計(jì)單元。根據(jù)原理構(gòu)建樣本序列[16]。具體過(guò)程為：首先計(jì)算不同窗口下單位面積上的地勢(shì)大小序列T（單位地勢(shì)梯度）（Ti＝平均起伏度（ti）/鄰域面積（si）檢驗(yàn)最有效，然后對(duì)序列T 取對(duì)數(shù)lnT，得到序列X（{xi，i＝2，3，4，…，23}））；利用公式（1）計(jì)算樣本序列X 的算數(shù)平均值，利用公式（2）和公式（3）分別計(jì)算樣本序列X 的統(tǒng)計(jì)量S 和Si的值。

1.2.1.3 劃分地形 根據(jù)數(shù)字地貌制圖規(guī)范，地形起伏度可分為7 級(jí)，即平原（＜30 m）、臺(tái)地（30～70 m）、丘陵（70～200 m）、小起伏山地（200～500 m）、中起伏山地（500～1 000 m）、大起伏山地（1 000～2 500 m）和極大起伏山地（＞2 500 m）[17]。山西地貌有平原、山地、臺(tái)地、丘陵、河漫灘五大一級(jí)類型。結(jié)合地形起伏度和地貌類型進(jìn)行地形分區(qū)。

1.2.2 Kriging 空間插值方法 Kriging 插值法是基于空間自相關(guān)性和二階平穩(wěn)假設(shè)，依據(jù)估算誤差最小和半方差函數(shù)分析，賦予一定鄰域內(nèi)樣點(diǎn)不同權(quán)重后，計(jì)算樣點(diǎn)值加權(quán)和即為插值結(jié)果。

式中，γ（h）為距離h 的樣點(diǎn)對(duì)的半方差；h 為2 個(gè)樣點(diǎn)之間的距離；N（h）表示樣點(diǎn)對(duì)距離為h 所有樣點(diǎn)對(duì)應(yīng)的個(gè)數(shù)；z（xi）和z（xi＋h）分別是采樣點(diǎn)有機(jī)質(zhì)在空間位置xi和xi＋h 處的實(shí)測(cè)值。當(dāng)γ（h）不與方向變化有關(guān)、只依賴于距離h 時(shí)，則稱γ（h）各向同性。

在數(shù)據(jù)量很大時(shí)需要在擬合變異函數(shù)之前進(jìn)行分組操作，這樣才可以任意距離下的變異函數(shù)進(jìn)行后續(xù)插值。球狀模型、指數(shù)模型、高斯模型、線性模型作為傳統(tǒng)理論擬合模型應(yīng)用比較廣泛，以下是對(duì)普通Kriging 插值原理進(jìn)行的說(shuō)明。

Kriging 插值預(yù)測(cè)結(jié)果的精確性是滿足無(wú)偏性和估計(jì)方差最小。

公式（5）和（6）經(jīng)推導(dǎo)，用變異函數(shù)來(lái)表示如下：

式中，k 為總的采樣點(diǎn)數(shù)；λi為克里金權(quán)重系數(shù)，表示每個(gè)采樣點(diǎn)xi處的有機(jī)質(zhì)含量對(duì)預(yù)測(cè)點(diǎn)x0的影響范圍；μ 為克里金拉格朗日乘子；γ（xi，xj）變異函數(shù)值是表示實(shí)際采樣點(diǎn)在空間距離xi和xj之間的變異值；γ（x0，xj）表示預(yù)測(cè)點(diǎn)x0與采樣點(diǎn)xj之間距離的變異值。

式中，z*（x0）表示預(yù)測(cè)點(diǎn)x0處的值；x1，…，xk為已知樣本的實(shí)際位置；z（x1），…，z（xk）為對(duì)應(yīng)的樣本有機(jī)質(zhì)實(shí)際含量。

1.2.3 預(yù)測(cè)結(jié)果準(zhǔn)確性檢驗(yàn) 本研究是通過(guò)回歸系數(shù)r 來(lái)驗(yàn)證采樣點(diǎn)位置處的預(yù)測(cè)精度，以及Isaaks 和Shrivastava 定義的平均預(yù)測(cè)誤差（ME）、平均標(biāo)準(zhǔn)誤差（ASE）、標(biāo)準(zhǔn)化均方根誤差（RMSSE）、均方根誤差（RMSE）[18]4 個(gè)驗(yàn)證指標(biāo)。其中，回歸系數(shù)r 值的大小是衡量變量之間的密切程度；均方根誤差RMSE 是觀測(cè)值與真實(shí)值偏差的平方和與觀測(cè)次數(shù)的比值的平方根來(lái)計(jì)算得到的。實(shí)際測(cè)量中，由于觀測(cè)樣點(diǎn)次數(shù)有限，所以用最接近的值來(lái)代替實(shí)際值。在一組數(shù)據(jù)中由于特殊值造成的誤差影響，會(huì)對(duì)標(biāo)準(zhǔn)誤差造成影響，所以預(yù)測(cè)精度用標(biāo)準(zhǔn)誤差來(lái)衡量。

2 結(jié)果與分析

2.1 鄰域分析與均值變點(diǎn)分析

利用Arcgis 對(duì)研究區(qū)在窗口大小分別為2，3，…，23 的情況下提取地形起伏度，窗口大小與平均起伏度的關(guān)系如表1 所示。

由表1 可知，平均起伏度值隨鄰域面積的增大逐漸增大，但當(dāng)鄰域面積達(dá)到一定的閾值時(shí)，其變化趨勢(shì)逐漸平緩。

平均起伏度和領(lǐng)域面積關(guān)系曲線符合邏輯斯蒂規(guī)律[16]。從圖2 可以看出，經(jīng)過(guò)統(tǒng)計(jì)學(xué)檢驗(yàn)，擬合度滿意，擬合方程為y＝20.171lnx＋14.056，擬合系數(shù)R2＝0.966 5，平均起伏度與鄰域面積成正比關(guān)系，趨勢(shì)由陡變緩處即為最佳統(tǒng)計(jì)單元面積。

S 和Si通過(guò)均值變點(diǎn)分析法計(jì)算得出，利用SPSS 統(tǒng)計(jì)軟件做出S 和Si差值的變化曲線如圖3所示。從圖3 可以看出，在第8 個(gè)點(diǎn)時(shí)S 和Si差值達(dá)到最大，該點(diǎn)是曲線由陡變緩的那一點(diǎn)，符合均值變點(diǎn)分析法的要求。由此通過(guò)鄰域分析法和均值變點(diǎn)分析計(jì)算得出，基于30 m 分辨率的dem 數(shù)據(jù)提取晉中市太谷區(qū)地形起伏度的最佳統(tǒng)計(jì)為9×9。

2.2 地形分區(qū)結(jié)果

本研究的最佳統(tǒng)計(jì)單元為9×9，在該窗口下該區(qū)域的地形起伏度范圍為0～241 m，按照數(shù)字地貌制圖規(guī)范[17]，該研究區(qū)耕地按地形起伏度可劃分為0～30 m 和30～70 m（圖4-A）。晉中市太谷區(qū)地貌類型為平原、臺(tái)地、丘陵、山地（圖4-B）。將地形起伏度分級(jí)結(jié)果與地貌類型進(jìn)行疊加分析，可分為0～30 m 的平原、臺(tái)地、丘陵、山區(qū)，30～70 m的平原、臺(tái)地、丘陵、山區(qū)共8 個(gè)區(qū)。

2.3 土壤有機(jī)質(zhì)統(tǒng)計(jì)分析

研究區(qū)土壤有機(jī)質(zhì)采樣點(diǎn)描述性統(tǒng)計(jì)如表2所示，利用SPSS 進(jìn)行分析檢驗(yàn)顯示，土壤有機(jī)質(zhì)含量符合正態(tài)分布規(guī)律。不同地形起伏度下地貌類型的土壤有機(jī)質(zhì)含量差別較大，符合Kriging 插值的條件；而且由于差異性也為接下來(lái)進(jìn)行半方差函數(shù)結(jié)構(gòu)建模以及依據(jù)地形分區(qū)插值提供依據(jù)。

表2 土壤有機(jī)質(zhì)含量描述性統(tǒng)計(jì)結(jié)果

由表3 可知，不同地形起伏度的地貌類型土壤有機(jī)質(zhì)含量的平均值、變異系數(shù)各不相同，從另一個(gè)角度證明了將地形起伏度和地貌類型作為輔助變量進(jìn)行分區(qū)Kriging 插值法的可行性。

表3 分區(qū)土壤有機(jī)質(zhì)含量描述性統(tǒng)計(jì)結(jié)果

2.4 地統(tǒng)計(jì)空間分析及Kriging 空間插值

利用地統(tǒng)計(jì)軟件GS+對(duì)各分區(qū)采樣點(diǎn)土壤有機(jī)質(zhì)含量進(jìn)行半方差函數(shù)及擬合參數(shù)計(jì)算[19-21]，其結(jié)果如表4 所示。

根據(jù)最優(yōu)擬合模型決定系數(shù)大、殘差小為依據(jù)，比較各模型參數(shù)，選擇最優(yōu)擬合指數(shù)參數(shù)[22]。根據(jù)塊金值的相關(guān)定義，原理上采樣點(diǎn)之間距離等于0 時(shí)，半變異函數(shù)值也是0，但由于存在測(cè)量誤差以及空間變異特征的影響，土壤樣點(diǎn)十分相近時(shí)，半變異函數(shù)值會(huì)發(fā)生變化，從而引起塊金值效應(yīng)，而且可以反映最小距離尺度下變量的變異特征及誤差。由表4 可知，不同地形起伏度和不同地貌類型下土壤有機(jī)質(zhì)與不分區(qū)有機(jī)質(zhì)的塊金值不同，說(shuō)明土壤有機(jī)質(zhì)的空間異質(zhì)性較大，分區(qū)與整體插值土壤有機(jī)質(zhì)空間變異不同。

表4 土壤有機(jī)質(zhì)含量的半方差函數(shù)模型和參數(shù)

由于該區(qū)域氣候差別很小，各分區(qū)的土壤有機(jī)質(zhì)結(jié)構(gòu)方差C 不同，所以該區(qū)域的空間變異特征主要是由地形因素造成的。C/（Co＋C）即空間相關(guān)度，是塊金值與基臺(tái)值的比值，該值反映了系統(tǒng)變量的空間相關(guān)性的大小。根據(jù)CAMBARDELLA 等于1994 年提出的劃分依據(jù)，塊金基臺(tái)比小于25%時(shí)，表明空間相關(guān)性強(qiáng)，大于75%則空間相關(guān)性較弱，處于二者之間，說(shuō)明具有中等強(qiáng)度的空間相關(guān)性。由表4 可知，在地形起伏度不同地貌類型土壤有機(jī)質(zhì)的C/（Co＋C）處于25%～75%時(shí)，表明土壤有機(jī)質(zhì)呈現(xiàn)中等強(qiáng)度的空間相關(guān)性。不同地形起伏度山地的C/（Co＋C）都為0.500，小于平原、臺(tái)地和丘陵地區(qū)，說(shuō)明山地上的耕地土壤有機(jī)質(zhì)受地形起伏度的影響更大。各分區(qū)的模型總體擬合效果較好，決定系數(shù)除地形起伏度為30～70 m 的山地之外都在0.72 以上，這與山地上的耕地少采樣點(diǎn)少有很大關(guān)系。地形分區(qū)下山地的變程相對(duì)較大，與晉中市太谷區(qū)山地地形地貌空間結(jié)構(gòu)相吻合。

與原始有機(jī)質(zhì)含量相比，不同地形起伏度下，平原、臺(tái)地、丘陵、山地的塊金值、基臺(tái)值都降低，丘陵的塊金基臺(tái)比升高，說(shuō)明不同分區(qū)土壤有機(jī)質(zhì)的變化比較小，相對(duì)比較均勻，表明該區(qū)域土壤有機(jī)質(zhì)含量進(jìn)行預(yù)測(cè)時(shí)考慮地形的合理性，把地形起伏度和地貌作為輔助變量對(duì)該研究區(qū)土壤有機(jī)質(zhì)含量預(yù)測(cè)的精度較好。

2.5 預(yù)測(cè)精度對(duì)比分析

由2 種方法所得到的實(shí)測(cè)值與預(yù)測(cè)值的散點(diǎn)分布圖和線性回歸方程如圖5 所示，以地形起伏度和地貌類型為輔助變量進(jìn)行分區(qū)的Kriging 插值得到的實(shí)測(cè)值與預(yù)測(cè)值的決定系數(shù)（0.463 4）明顯高于普通的Kriging 插值法得到的（0.261 5）。

ME 越接近0 越滿足預(yù)測(cè)的無(wú)偏性；RMSSE 越接近1 說(shuō)明預(yù)測(cè)標(biāo)準(zhǔn)誤差越準(zhǔn)確；RMSE 和ASE 盡可能小，預(yù)測(cè)值與實(shí)測(cè)值的差值更小，使結(jié)果更精確[23-25]。從表5 可以看出，地形起伏度在0～30 m 山地的無(wú)偏性最好；在30～70 m 的丘陵和山地預(yù)測(cè)值與實(shí)測(cè)值的偏差最小，所以結(jié)果更準(zhǔn)確。各分區(qū)土壤有機(jī)質(zhì)的交叉驗(yàn)證均滿足上述條件[26]。

2.6 土壤有機(jī)質(zhì)的空間分布特征

通過(guò)普通Kriging 和分區(qū)Kriging 這2 種方法預(yù)測(cè)得到的土壤有機(jī)質(zhì)空間分布如圖6 所示，分區(qū)Kriging 插值結(jié)果空間分布特征及空間遞變規(guī)律更加明顯，普通Kriging 插值得到的土壤有機(jī)質(zhì)含量范圍為11.47～22.41 g/kg，而分區(qū)Kriging 插值得到的土壤有機(jī)質(zhì)含量為7.38～25.95 g/kg?？傮w上看，研究區(qū)2 種方法估計(jì)的土壤有機(jī)質(zhì)空間分布格局基本一致，都表現(xiàn)為西北部的土壤有機(jī)質(zhì)含量較高，東南部較低。

從2 種預(yù)測(cè)結(jié)果局部分布特征看，土壤有機(jī)質(zhì)含量差異還是比較明顯，分析發(fā)現(xiàn)，普通Kriging 插值，由于平滑效應(yīng)使不同地形起伏度地貌類型的有機(jī)質(zhì)含量差異變低[27-28]，插值結(jié)果不能精確反映出局部信息，例如，山地、丘陵上耕地由于個(gè)別采樣點(diǎn)的土壤有機(jī)質(zhì)含量偏高，造成該區(qū)域的插值結(jié)果偏高，而且分布規(guī)律連續(xù)規(guī)整，因此無(wú)法體現(xiàn)出太谷區(qū)的局部信息；并且普通Kriging 插值并沒(méi)有對(duì)采樣點(diǎn)根據(jù)地形歸類，相鄰距離近的點(diǎn)忽略了空間位置因素造成的誤差，且普通插值法忽略了土壤有機(jī)質(zhì)空間平穩(wěn)過(guò)渡，造成預(yù)測(cè)的土壤有機(jī)質(zhì)分布特征不明顯。但分區(qū)Kriging 是按照地形分區(qū)對(duì)采樣點(diǎn)進(jìn)行歸類，可以有效消除鄰近點(diǎn)插值誤差，更能反應(yīng)土壤有機(jī)質(zhì)的空間特點(diǎn)以及研究區(qū)域局部分布情況，更適合用于分析不同因素與土壤有機(jī)質(zhì)空間規(guī)律之間的關(guān)系。

3 結(jié)論

本研究結(jié)果表明，結(jié)合地形起伏度和地貌類型進(jìn)行分區(qū)Kriging 插值法得到的研究區(qū)土壤有機(jī)質(zhì)含量的范圍為7.38～25.95 g/kg，而直接Kriging 插值法得到的研究區(qū)土壤有機(jī)質(zhì)含量范圍為11.47～22.41 g/kg，2 種方法所得范圍大致一致，但分區(qū)Kriging 插值法預(yù)測(cè)精度更高。這是因?yàn)檠芯繀^(qū)地貌和地形起伏度不同，導(dǎo)致區(qū)域內(nèi)土壤有機(jī)質(zhì)含量不同，結(jié)合地形起伏度和地貌進(jìn)行分區(qū)以及對(duì)采樣點(diǎn)進(jìn)行歸類可以消除鄰近采樣點(diǎn)對(duì)其空間插值預(yù)測(cè)精度的影響，從而提高其預(yù)測(cè)精度。基于地形分區(qū)的Kriging 插值法對(duì)丘陵和山地的預(yù)測(cè)結(jié)果更好，因此，對(duì)于研究地形復(fù)雜山區(qū)的有機(jī)質(zhì)空間變異規(guī)律，可以考慮借助輔助變量對(duì)空間規(guī)律進(jìn)行精細(xì)研究。