基于阻隔因素的耕地質量分等因素插值方法研究

楊永俠 郭雅萍 張 函 張麗紅 桑 婧

(1.中國農業(yè)大學土地科學與技術學院， 北京 100083; 2.自然資源部農用地質量與監(jiān)控重點實驗室， 北京 100035)

0 引言

耕地質量調查評價是新時期下對耕地數量、質量和生態(tài)“三位一體”保護的重要基礎，同時也是有效保護耕地資源、保障國家糧食安全、維持農業(yè)可持續(xù)發(fā)展的重要基石，可為解決國家耕地的占補平衡、基本農田劃定、土地利用規(guī)劃等問題提供科學依據[1-3]。

耕地質量評價是以縣為單位劃分耕地質量分等單元，縣域內動輒有過萬分等單元，其數量大、分布廣，所以分等因素一般通過樣點(樣地)數據采用空間插值的方法獲得。對于平原地區(qū)或者分等因素變化不大的區(qū)域，一般插值方法可以較好地擬合實際情況，但當區(qū)域內地形復雜、海拔落差較大時，利用克里金或者反距離權重插值的方法不能很好地擬合實際情況。目前，對于此類情況主要是通過增加調查樣點提高插值的準確性，但是位于山地丘陵區(qū)耕地的地形變化過于復雜，增加控制樣點的方法同樣也會消耗過多的人力物力，且插值精度的提高并不明顯。

在復雜地區(qū)土壤屬性空間預測方面，國內學者對空間預測模型精度的提高[4-6]及預測性制圖[7-12]等進行了研究，并已取得多方面的研究進展。經研究證實，利用輔助變量可以提高土壤屬性的空間預測精度，但是不能很好再現原始區(qū)域變量的空間結構及變量之間的相互關系，而且不能進行定量的不確定性分析，這使與輔助變量結合的預測結果應用受到一定的限制。

基于以上問題，本文在耕地質量評價體系的基礎上，提出山區(qū)分等因素基于阻隔因素的插值估算方法，以期提高耕地質量等級評價結果的準確性，為第三次全國土地調查耕地等別調查評價工作提供借鑒。

1 研究區(qū)和數據

1.1 研究區(qū)概況

研究區(qū)選擇青海省海東市平安區(qū)湟水流域中游南側，海東市中心腹地，西距西寧市35 km處，地理坐標為東經101°49′～102°10′、北緯36°15′～36°34′之間，境內南北長約33.3 km，東西寬約23 km，總面積769 km2。

平安區(qū)屬大陸性氣候，春季干旱多風，夏季涼爽，冬季寒冷。年平均氣溫7.6℃，年內降雨分布不均，無霜期218 d。平安區(qū)境內多為山區(qū)，且有東溝、六道溝、老虎溝等12條溝岔，成為山巒起伏、山川相間的主要河谷地帶。平安區(qū)境內的大部分地區(qū)海拔在2 038～4 127 m之間(圖1)，起伏連綿的地勢以及復雜多樣的氣候，構成了農業(yè)生產在空間上比較明顯的地區(qū)差異，最終形成了川水、淺山、腦山3種地貌類型區(qū)。

圖1 平安區(qū)海拔分布Fig.1 Elevation map of Ping’an District

1.2 數據來源與處理

采用數據來源于青海省平安區(qū)2007—2011年5年的測土配方施肥項目的調查數據，在綜合分析本研究區(qū)內耕地分布特點的基礎上，均勻選取了280個調查樣點數據，利用ArcGIS中子集要素工具(Subset Features)隨機確定240個點作為插值點，40個點作為驗證點進行精度評估(圖2)。

圖2 插值樣點分布Fig.2 Interpolation sample distribution map

2 研究方法

2.1 可插值分等因素分析

《農用地質量分等規(guī)程》[13]中指出，農用地質量分等因素分為推薦因素和自選因素兩類。推薦因素由國家統(tǒng)一確定，并且分區(qū)、分地貌類型給出，自選因素由省級土地行政主管部門確定[14]。總的來說推薦因素和自選因素可以分為4種類型，即土壤條件、地貌條件、水文類以及農田基本建設類，每種類型所涉及的分等因素如表1所示。

平安區(qū)在農用地質量分等過程中選用的分等因素有：表層土壤質地、地形坡度、灌溉保證率、土壤pH值、土壤有機質含量、有效土層厚度。由于地形坡度這個分等因素的空間相關性較差且不滿足地理學第一定律；而表層土壤質地和灌溉保證率的屬性值在空間上也不連續(xù)，故均不適合利用插值對未知樣點值進行估算；因此本文考慮到分等因素的含義、屬性值類型、插值方法，選取了平安區(qū)有機質含量、土壤pH值和有效土層厚度3個分等因素進行插值方法研究。雖僅根據3個因子不能對耕地質量進行分等定級，但這3個因子在山地丘陵區(qū)耕地質量評價時所占的權重均比較大，且其屬性值會很大程度影響最終的評價結果，故需保證這三者的準確性。

表1 分等因素統(tǒng)計Tab.1 Statistics of factors

2.2 傳統(tǒng)的空間插值方法

2.2.1反距離權重插值

反距離權重法是空間插值的常用而又簡單的方法。它以樣本點間距離為權重進行加權平均，距離插值點越近的樣本所占權重越大，這種方法簡單方便，效果直觀而且效率較高[15]。但如果樣本中有極端值則會影響較大范圍的插值結果。其計算公式為

(1)

反距離權重法和克里金法的計算公式相同，只是在反距離權重法中權重λi僅取決于預測位置的距離。

2.2.2樣條函數插值

樣條函數法是指通過使用最小化整體表面曲率的函數來估計值，生成恰好經過樣本點的平滑表面，如同將一個軟膜插入并經過各個已知樣點，并且表面的總曲率最小[17]。此方法最適合生成平緩變化的表面，比如海拔、地下水位等。其計算公式為

(2)

式中S——所求點函數值

T——原始點函數值

R——距離函數N——樣點數

λj——通過求解線性方程組而獲得的系數

rj——點(x,y)與點j之間的距離

2.3 基于阻隔因素的插值方法

在山地丘陵區(qū)利用插值對未知樣點進行估值時，區(qū)域內極有可能會有較大的山體或者河流分布，對插值效果影響較大，所以在進行插值分析時應考慮這一因素。已有的研究表明，土壤屬性具有高程地帶性的分布規(guī)律，并且可以用高程作為輔助變量對土壤的屬性進行模擬預測，但是此類研究側重于對原數據的優(yōu)化，不能反映數據本身的結構與特征；因此本文在不改變原數據的狀態(tài)下，從模擬預測過程中引入阻隔因素，便于更好地貼合土壤屬性的實際分布狀況。

阻隔因素是指在研究區(qū)中會中斷表面連續(xù)性的線狀要素的位置。山地丘陵區(qū)實際工作中阻隔因素可能為山脊、山谷、河流，這樣的地形因素存在時間長，對周圍土地影響大，會使阻隔因素兩側的數據不連續(xù)，擬合的數據曲面容易發(fā)生突變。尤其山體阻隔影響，迎風坡與背風坡風化差異，陽坡與陰坡光照輻射差異，氣流和能量對土壤變化影響尤為顯著。而且在山地丘陵區(qū)，耕地多分布于山谷和河谷兩側，所以本研究選取海拔落差較大的山體與大型河流作為阻隔因素進行研究。

基于阻隔因素的插值方法具體步驟為：首先提取阻隔因素，本文中阻隔因素提取涉及河流和山體。大型河流分布數據可以通過實地調查或者河流分布圖獲得，方法較為簡單，本文不做贅述；山體阻隔因素提取，山脊線的獲取采用水文分析的方法[16]，算法思路為：山脊線同時也是分水線，對于在分水線上的柵格，是水流的起源，通過地表徑流模擬計算，應當只有流出方向的表面徑流，沒有流入方向，也就是此處柵格的匯流量為零。提取匯流值為零的柵格，就可以得到山脊線。然后提取不穿越耕地，而且其兩側距離耕地較遠的山脊線為山體阻隔線。通過矢量化、連接、平滑等操作確定分水線的位置。利用ArcGIS軟件對DEM計算地表徑流，提取匯流量為零的柵格，經過篩選、矢量化、平滑得到平安區(qū)阻隔因素分布圖，如圖3所示，將提取的阻隔因素加入插值的過程中。ArcGIS軟件中的反距離權重插值工具提供“輸入折線障礙要素”選項，此處的折線要素是指在搜索輸入采樣點時用作中斷或限制的折線要素，因此可以將提取的阻隔因素作為其折線要素，其他參數采用默認設置進行插值分析；ArcGIS軟件也提供了“含障礙的樣條函數”插值工具，同樣將阻隔因素作為障礙要素進行輸入，其他參數保持不變。

圖3 阻隔線分布圖Fig.3 Block line distribution map

2.4 精度檢驗方法

采用交叉驗證[18-20]對反距離權重插值、樣條函數插值、基于阻隔因素的反距離權重插值、基于阻隔因素樣條函數插值這4種插值方法的精度進行對比分析，采用平均絕對誤差(MAE)、平均相對誤差(MRE)、均方根誤差(RMSE)、相對均方根誤差(RMSEr)4個評判標準來評價插值結果。其中，平均絕對誤差能夠估算出獲取的估計值的誤差范圍；平均相對誤差可以判斷出估計值和實測值之間的誤差；均方根誤差反映利用樣點數據的估值靈敏度和極值效應，其值越小越好；相對均方根誤差越小則誤差越小，精度越高[21-23]，計算公式如下

(3)

(4)

(5)

(6)

式中xi1——驗證樣點的測量真值

xi2——驗證樣點預測值

n——驗證樣本總數

3 結果與分析

3.1 原始數據描述性統(tǒng)計分析

對原始樣本的有效土層厚度、有機質含量、土壤pH值進行基本的描述性統(tǒng)計，從表2可知，有效土層厚度的范圍為27～153 cm，且極差最大；有機質質量比范圍是9～62 g/kg，土壤pH值極差最小；偏度和峰度是描述數據是否符合正態(tài)分布的2個重要指標，當二者都為0時，表明數據符合標準的正態(tài)分布。當偏度大于0時，表明數據集中分布在左側，向右延伸，小于0時，集中分布在右側，向左延伸。當峰度大于0時，整個數據分布形態(tài)比標準的正態(tài)分布高聳，數據多集中分布在平均值附近，小于0時，整個數據分布形態(tài)比標準正態(tài)分布平坦，數據分布較為分散。由表2可知，3個分等因素均不符合標準的正態(tài)分布；變異系數是反映數據分布狀態(tài)的一項重要指標，主要反映數據的離散程度。由表2可知，有效土層厚度和有機質含量屬于中等變異性，pH值屬于弱變異性，所以有效土層厚度和有機質含量較pH值更易受環(huán)境變量的影響而變化。

表2 原始樣本描述性統(tǒng)計特征值Tab.2 Descriptive statistical eigenvalues of original sample

3.2 插值結果分析

3.2.1土壤有機質含量結果分析

調查樣點有機質含量通過反距離權重、基于阻隔因素的反距離權重、樣條函數、基于阻隔因素的樣條函數4種插值方法獲得的插值結果如圖4所示。

圖4 有機質含量插值結果對比Fig.4 Comparisons of organic matter content interpolation results

土壤有機質不僅是植物養(yǎng)分供給的源泉之一，而且是保持土壤良好的物理性質的物質。平安區(qū)全區(qū)耕地土壤的有機質質量比的范圍在9～62 g/kg之間，就全區(qū)而言，平安區(qū)的有機質水平均較低。由圖4可以看出，這4種不同的插值方式均能較好地模擬平安區(qū)土壤有機質的分布狀況，呈現出由北向南逐漸遞增的趨勢。

圖5 土壤pH值插值結果對比Fig.5 Comparison of soil pH value interpolation results

4種插值方法得到的有機質含量的空間分布相一致，但其插值結果卻存在一定的差異，通過反距離插值、樣條函數插值、基于阻隔因素的反距離權重、基于阻隔因素的樣條函數插值得到的有機質含量范圍分別是9.065 94～61.966 1 g/kg、8.989 31～73.701 6 g/kg、9.064 9～61.966 1 g/kg、8.896 7～72.529 7 g/kg，通過對比可知，反距離權重插值的結果與有機質真實值的范圍最為接近；從平滑度和局部變異性來看，反距離權重插值出現了“牛眼”現象，樣條函數插值結果圖較平滑、連續(xù)；且對比基于阻隔因素與未基于阻隔因素的反距離插值和樣條函數插值結果，可以明顯看出，基于阻隔因素后對其兩側數據的插值產生的影響，能更真實反映土壤有機質分布情況的細節(jié)特征，因此基于阻隔因素的插值方法適用性更強。

對40個實際調查樣點數據的插值精度采用交叉驗證的方法進行分析，從表3可看出，基于阻隔因素的樣條函數插值方法精度最高，比樣條函數插值的平均絕對誤差、平均相對誤差、均方根誤差、相對均方根誤差分別提高了11.23%、10.98%、7.54%和9.20%；比反距離權重插值方法的平均絕對誤差、平均相對誤差、均方根誤差、相對均方根誤差分別提高了31.27%、25.60%、31.74%和21.39%。反距離權重、樣條函數、基于阻隔因素的反距離權重和基于阻隔因素的樣條函數4種插值方式的實測值與預測值的決定系數分別為：0.876、0.920、0.878和0.927；綜上可知基于阻隔因素的樣條函數插值方式最優(yōu)。

表3 有機質含量插值結果精度Tab.3 Accuracy analysis of organic matter content interpolation results

3.2.2土壤pH值結果分析

調查樣點土壤pH值通過反距離權重、基于阻隔因素的反距離權重、樣條函數、基于阻隔因素的樣條函數4種插值方法獲得的插值結果，如圖5所示。

平安區(qū)的全區(qū)耕地的土壤pH值范圍在7.83～8.56之間，呈現出從中部分別向南北遞減的分布趨勢。由圖5可看出，4種插值方法均能較好地反映平安區(qū)全區(qū)的土壤pH值的空間分布特征。

4種插值方法得到的土壤pH值的空間分布相一致，但其插值結果范圍卻存在一定的差異，通過反距離插值、樣條函數插值、基于阻隔因素的反距離權重、基于阻隔因素的樣條函數插值得到的土壤pH值范圍分別是：7.831 34～8.539 96、7.256 99～8.645 68、7.830 00～8.539 96、7.758 4～8.553 67，通過對比可知，反距離權重法對pH值估算的結果與真實值的范圍最為接近；但是采用反距離權重進行插值的結果在插值點高值附近會出現“牛眼”，這是由于反距離權重插值是一種精確插值方式，不改變插值點原始值；耕地質量調查評價工作要求，在插值過程中應該保證調查樣點的值不變，才能保證耕地質量評價的準確性；分別對比基于阻隔因素與未添加阻隔因素的反距離插值和樣條函數插值結果，可以得知基于阻隔因素的插值結果反映細節(jié)比較清晰，因此基于阻隔因素的反距離插值方法對土壤pH值進行插值的適用性最優(yōu)。

從表4可知，基于阻隔因素的反距離權重插值方法精度最高，比反距離權重插值的平均絕對誤差、平均相對誤差、均方根誤差、相對均方根誤差分別提高了2.94%、2.56%、1.85%和1.90%；比樣條函數插值方法的平均絕對誤差、平均相對誤差、均方根誤差、相對均方根誤差分別提高了11.19%、11.63%、14.39%和14.88%；反距離權重、樣條函數、基于阻隔因素的反距離權重和基于阻隔因素的樣條函數4種插值方式的實測值與預測值的決定系數分別為0.527、0.560、0.554和0.476；雖然樣條函數和基于阻隔因素的反距離插值方法兩者的決定系數均比較高，且相差不大，但是后者插值精度更高，故研究區(qū)內對土壤pH值最適合的插值方法是基于阻隔因素的反距離權重插值方法。

表4 土壤pH值插值精度Tab.4 Soil pH value interpolation accuracy

3.2.3有效土層厚度結果分析

調查樣點有效土層厚度通過反距離權重、基于阻隔因素的反距離權重、樣條函數、基于阻隔因素的樣條函數4種插值方法獲得的插值結果如圖6所示。

圖6 有效土層厚度插值結果對比Fig.6 Comparison of effective soil thickness interpolation results

平安區(qū)耕地的有效土層厚度在27～153 cm之間，其分布特征是：東北部和西南部最低，西北部的有效土層厚度最高；且呈現出從中部向東西兩側遞減的趨勢。由圖6可以看出，4種插值方法均能較好地反映土壤有效土層厚度的空間分布特征。

4種插值方法得到的有效土層厚度的空間分布相一致，但其插值結果卻存在一定的差異，通過反距離權重、樣條函數插值、基于阻隔因素的反距離權重、基于阻隔因素的樣條函數插值得到的有效土層厚度范圍分別是：27.567～152.824 cm、26.002 7～204.752 0 cm、27.567～152.872 cm、26.147 3～166.074 0 cm，通過對比可知，反距離權重結果與真實值的范圍最為接近，樣條函數插值次之；同樣采用反距離權重進行插值的結果與土壤pH值類似，在插值點高值附近會出現“牛眼”，且其空間分布圖過渡平滑性欠佳，可觀性較差；樣條函數插值后的結果較好，圖形表面平滑漸進；且分別對比基于阻隔因素與未基于阻隔因素的反距離權重和樣條函數插值結果圖，基于阻隔因素的插值結果反映了阻隔因素對于其兩側影響的細節(jié)，形成明顯的分界線，其插值總體保留了豐富的細節(jié)，還反映了分等因素的空間漸變特征，更符合其實際的狀況，因此基于阻隔因素的樣條函數插值方法對土壤有效土層厚度的適用性最優(yōu)。

從表5可知，基于阻隔因素的樣條函數插值方法精度最高，比樣條函數插值的平均絕對誤差、平均相對誤差、均方根誤差、相對均方根誤差分別提高了15.08%、11.74%、17.41%和9.40%；比反距離權重插值方法的平均絕對誤差、平均相對誤差、均方根誤差、相對均方根誤差分別提高了26.33%、23.03%、25.05%和24.63%；反距離權重、樣條函數、基于阻隔因素的反距離和基于阻隔因素的樣條函數4種插值方式的實測值與預測值的決定系數分別為0.869、 0.852 、0.865和0.901；故基于阻隔因素的樣條函數的插值方法最適宜有效土層厚度的區(qū)域空間預測。

表5 有效土層厚度插值精度Tab.5 Interpolation accuracy of effective soil thickness values

4 討論

4.1 基于阻隔因素插值結果與耕地質量年度更新評價結果對比

將土壤有機質含量、土壤pH值和有效土層厚度3個分等因素插值結果與耕地質量評價成果中因素指標值的值域和均值對比可得，分等因素的值域與均值均發(fā)生明顯變化(表6)。基于阻隔因素的樣條函數插值得到的耕地有機質含量為11.84～35.61 g/kg，均值為19.63 g/kg，比年度評價成果均值26.45 g/kg低6.82 g/kg，其平均絕對誤差、平均相對誤差、均方根誤差及相對均方根誤差分別提高了 87.35%、87.69%、87.78%和88.52%；基于阻隔因素的樣條函數插值得到的耕地有效土層厚度為43.74～149.15 cm，均值為73.59 cm，比年度評價成果均值76.25 cm低2.66 cm，其平均絕對誤差、平均相對誤差、均方根誤差及相對均方根誤差分別提高了51.88%、52.80%、63.18%和61.89%；基于阻隔因素的反距離權重插值得到土壤pH值的值域為7.98～8.37，均值為8.26，比年度評價成果均值8.28低0.02，其平均絕對誤差、平均相對誤差、均方根誤差及相對均方根誤差分別提高了72.32%、72.36%、69.19%和69.25%。3種分等因素利用基于阻隔因素得到的結果，耕地圖斑地塊的分等因素值均有減少，但總體來看，值域和均值變化不大，且對土壤有機質、有效土層厚度和土壤pH值3個分等因素的插值精度均有較大提高。在本文研究中發(fā)現，基于阻隔因素的確定性插值方法對不同的土壤性質表現不一致，土壤有機質和有效土層厚度最適宜的插值方式是基于阻隔因素的樣條函數插值，而在土壤pH值插值中，基于阻隔因素的反距離權重插值精度最高，這是由于反距離權重和樣條函數插值原理不同和3個分等因素原始數據的特點各異造成的。由表2可知，土壤pH值屬于弱變異性，方差最小，相較于另外兩個分等因素其數據分布較為集中；而反距離權重插值由于其根據插值點與樣本點之間的距離為權重進行加權平均，在樣本分布較為集中、均勻時效果最佳。

表6 插值結果與年度更新數據對比Tab.6 Comparison of interpolation results and annual update data

4.2 插值結果的進一步修正

在實際情況中對于山體、河流而言，迎風坡與背風坡風化存在著差異，陽坡與陰坡光照的輻射也不相同，氣流和能量對土壤變化影響也尤為顯著。而且在山地丘陵區(qū)，耕地多分布于山谷和河谷兩側。土壤有機質、土壤pH值和有效土層厚度的實際分布特征不僅會受到地形地貌的影響，也會受到其他自然因素與人為因素的影響。土壤有機質的具體分布會受到高程[24]、坡度、地形濕度指數、土壤類型[25]等自然因素和種植模式[26]、施肥結構[27]、灌溉、土地利用方式等人為因素的影響，其中高程在大尺度上與土壤有機質的關系最為明顯，土壤類型主要是由于成母土質的差異影響土壤有機質的分布[28]。在河流流經的地區(qū)土壤質地對土壤有機質的空間分布也有顯著的影響[29]；土壤的pH值影響土壤中某些養(yǎng)分元素的有效性，從而影響植物生產，而施用肥料也會對土壤pH值產生影響，相關研究表明土壤pH值與同樣也會受到地形地貌、成母土質[30]、土壤類型[31]等自然因素和農耕、施肥措施、工業(yè)發(fā)展等人為因素的影響；土層厚度關系到植物的扎根條件，而且與養(yǎng)分和水分的保蓄能力相關，薄層土壤對深根性植物的生長有限制，而有效土層厚度和地形地貌、土壤類型、土地利用類型、巖石裸露率、植被覆蓋度有一定的相關性。

對平安區(qū)有機質含量、有效土層厚度以及土壤pH值3個分等因素的屬性值采用基于阻隔因素的方法對其進行了插值對比。但是分等因素種類繁多且分布類型不一，各個地區(qū)采用的分等因素也不盡相同，而本文僅對山地丘陵區(qū)耕地質量分等因素屬性值的預測提出了新的方法，并非能適用于所有地區(qū)的所有分等因素；而且本文研究的范圍僅是小尺度(縣區(qū)級)，若獲取更為準確的大尺度范圍內(地形復雜區(qū))的分等因素屬性值，還需要在本方法的基礎上考慮多重環(huán)境因子。

5 結論

(1)采用不同的插值方法對土壤有機質含量、土壤pH值、有效土層厚度3種分等因素進行插值后，3種分等因素的插值結果的空間分布趨勢的模擬與其實際分布基本一致。

(2)由基于阻隔因素的插值方法與不添加阻隔因素的插值方法進行對比可知，基于阻隔因素的插值結果反映了阻隔因素對于其兩側影響的細節(jié)，形成明顯的分界線，其插值總體保留了豐富的細節(jié)，還反映出分等因素的空間漸變特征，更符合其實際的狀況，同時也具有更高的預測精度。

(3)對研究區(qū)的土壤pH值進行插值，基于阻隔因素的反距離權重插值方法優(yōu)于其他3種插值方法。

(4)對研究區(qū)的土壤有機質含量和有效土層厚度進行插值，基于阻隔因素的樣條函數插值方法的精度均明顯優(yōu)于其他3種插值方法，且二者實測值與預測值之間的決定系數分別為0.927和0.901；對土壤有機質含量的插值，基于阻隔因素的樣條函數插值方法比樣條函數插值的平均絕對誤差、平均相對誤差、均方根誤差、相對均方根誤差分別提高了11.23%、10.98%、7.54%和9.20%；對有效土層厚度的插值，基于阻隔因素的樣條函數插值方法比樣條函數插值的平均絕對誤差、平均相對誤差、均方根誤差、相對均方根誤差分別提高了15.08%、11.74%、17.41%和9.40%。