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

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

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

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

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

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

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

1 材料和方法

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

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

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

1.2 研究方法

1.2.1 地形分區(qū)

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

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

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

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

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

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

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

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

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

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

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

2 結(jié)果與分析

2.1 鄰域分析與均值變點分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.5 預(yù)測精度對比分析

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

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

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

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

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

3 結(jié)論

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