表征土壤有機碳區(qū)域分布的優(yōu)化空間插值模型研究——以福建省龍海市為例

范勝龍，黃炎和，林金石

（福建農林大學 資源與環(huán)境學院，福州 350002）

土壤有機碳是土壤理化性質中最重要、最具代表性的因子，其在空間分布上具有不均一性，是發(fā)生變化的連續(xù)體［1－2］。對土壤有機碳儲量、分布、轉化、衰減機理進行研究，并揭示其影響因素和生態(tài)效應是目前土壤有機碳研究的熱點問題。而土壤樣點的代表性是研究結果是否準確的前提，大規(guī)模的土壤屬性空間信息的采集需要花費大量的人力、物力和財力。因此，研究如何利用有限的樣本數據來獲得更為詳盡的土壤屬性空間分布信息的方法具有重要意義。

空間插值模型是實現土壤有機碳含量從離散的點狀信息向面狀連續(xù)信息轉換的有力工具，是表征土壤有機碳空間分布特征的重要手段。目前，在土壤屬性空間預測的研究方面，地統計學被證明為最有效的空間插值方法［3－6］。綜觀學者們的研究，克里格插值是使用最多的空間預測方法，但在地形復雜地區(qū)，普通克里格方法的應用有一定的局限性。輔助信息用以提高目標變量的預測精度，在土壤學的研究中已經達成共識。協同克里格法就是一種利用輔助的方法，并在土壤科學中得到了廣泛的應用［7］。但相對于縣級或更大研究區(qū)域而言，為了達到優(yōu)化插值目的，而多調查另一數據，無疑增加了調查成本，且通過調查所得兩組數據一般樣點數目相同，因此，采用協同克里格法增加了樣點調查成本。本研究以福建省龍海市為典型區(qū)，以耕地中耕層土壤有機碳為研究對象，設計結合地貌類型、土壤類型和土地利用類型等信息的克里格插值模型，研究縣級尺度上土壤有機碳空間預測的優(yōu)化空間插值模型及其與樣點密度的關系。

1 材料與方法

1.1 研究區(qū)概況

龍海市位于福建省漳州市東部沿海，全市總面積1 289.72 km2，屬南方丘陵區(qū)。境內地勢北部、西部、南部三面環(huán)山、中部為平原、東南部臨海。土壤形成受成土母質的影響極大，由于母質分布的區(qū)域性不同，導致不同區(qū)域土壤變化過程具有明顯的差異性，加上長期以來自然因素和人為因素的作用，形成土壤類型的多樣性。根據第二次土壤普查資料，龍海市土壤分為6個土類，16個亞類，52個土屬，70個土種。

1.2 研究方法

1.2.1 土壤樣品采集及測定方法 考慮到格網法是有關土壤屬性調查樣點布設的最常用方法且能夠較為準確地表征土壤屬性的空間變異［8－10］，本研究采用格網法布設樣點。根據前人相關研究成果和研究區(qū)域的大小，本研究采用0.5 km×0.5 km（X1）、1 km×1 km（X2）、2 km×2 km（X3）、4 km×4 km（X4）四種格網布設樣點。

土壤樣品采集的具體做法是在樣點附近20 m的范圍內取5個耕作層（0—20 cm）土樣，然后混合成一個土壤樣品，用四分法取1.0 kg帶回實驗室，采用常規(guī)的低溫外加熱重鉻酸鉀氧化－滴定法測定土壤有機碳含量。在采樣的同時，用GPS記錄每個樣點的經緯度信息，并描述各樣點的土壤、土地利用類型、地貌及相關環(huán)境信息。所有土壤樣品采集在2009年11—12月農作物收割完成后進行。

1.2.2 土壤有機碳空間插值模型 由于地貌類型、土壤類型和土地利用類型對土壤有機碳空間分布均有重要影響，可以將其作為提高土壤屬性預測精度的輔助信息［11］。為研究上述輔助信息對優(yōu)化空間插值結果的影響，本研究在采用格網法取樣的基礎上，分別設計了結合地貌類型信息的克里格方法（KDM）、結合土地利用類型信息的克里格方法（KDL）、結合土壤類型信息的克里格方法（KTR）、結合地貌—土壤類型信息的克里格方法（KDMTR）、結合土地利用類型—土壤類型信息的克里格方法（KDLTR）五種結合類型信息的方法，并與直按采用各樣點土壤有機碳含量數據進行普通克里格插值（KYJZ）的結果進行比較分析。上述結合類型信息的方法將每一個樣點的土壤有機碳含量值z（xkj）分為相同類型均值μ（tk）和殘差r（xkj）之和。用公式表示為：

式中：z（xkj）——樣品的土壤有機碳含量；μ（tk）——相同類型樣品的均值；r（xkj）——樣品土壤有機碳含量與其相同類型樣品的均值之差，稱為“殘差”。

研究將殘差作為一個新的區(qū)域變量r（xkj）進行普通克里格插值，空間插值利用ArcGIS軟件中的地統計分析模塊完成。各待估點的土壤有機碳含量預測值z＊（xkj）為類型均值μ（tk）與殘差估計值r＊（xkj）之和。

本研究中，TR法、DMTR法和DLTR法中的土壤類型均劃分到土屬級別。土壤類型資料來自于第二次土壤普查成果中的1∶5萬土壤圖和《龍海土壤》；土地利用現狀資料來源于第二次全國土地調查標準時點統一更新數據庫；地貌劃分來源于福建省農用地分等成果更新項目成果中的指標區(qū)圖。

1.2.3 預測結果檢驗 為了檢驗各種空間插值方法的效果，本研究在全市范圍內隨機選擇259個驗證樣點，驗證樣點包含了預測樣點中所有地類（土壤類型和地貌類型），且每種類型的驗證樣點數量不少于3個。同時，在選擇驗證樣點時也兼顧空間分布的均勻性。通過驗證點的實測值與該驗證點通過各種插值方法所得的預測值的相關系數（r）及其均方根誤差（RMSE）來評價其預測精度的高低。r值越大、RMSE越小則預測精度越高，反之精度越低。

1.2.4 樣點密度對土壤有機碳空間插值模型精度的影響研究 在選定空間插值方法的基礎上，為探討樣點密度與土壤有機碳空間預測精度之間的量化關系，本研究利用布設的4種不同密度格網所得樣點，采用設計的6種結合類型信息的克里格法分別進行空間預測，并通過各驗證點實測值與預測值的RMSE，對其預測精度進行驗證，研究滿足精度要求的空間插值模型及其所需樣點密度。

2 結果與分析

2.1 龍海市土壤有機碳的空間分布特征

本研究采用4種不同密度格網疊套于龍海市1∶1萬土地利用現狀圖中進行樣點布設，全市共獲得1 133個耕地圖斑上的預測樣點，在所有樣點中，土壤有機碳含量最小值為3.50 g／kg，最大值為86.70 g／kg，兩者相差了近23倍，土壤有機碳含量均值為24.45 g／kg，變異系數43.35%，說明南方丘陵區(qū)縣級尺度上農田土壤有機 碳含量變異較大。

2.2 結合不同類型信息的土壤有機碳含量空間插值

利用ArcGIS軟件對4種格網密度下耕層有機碳含量數據去除各種類型均值后殘差采用普通克里格法進行空間插值，結果如圖1所示。受篇幅的限制，本研究僅以0.5 km×0.5 km格網為例，列出結果圖并進行分析。

圖1 不同方法所得插值結果

從圖1中可以看出，各種方法所得插值圖反映的土壤有機碳的空間分布趨勢基本一致，即北部和中部九龍江下游沖積平原的土壤有機碳含量較高，而東部港尾和隆教兩鄉(xiāng)鎮(zhèn)沿海的太武山區(qū)土壤有機碳含量最低。各種方法所得插值圖反映的土壤有機碳的這種分布趨勢與龍海市土壤有機碳實際的空間分布情況基本一致，表明所用各種空間插值方法能基本準確反映出研究區(qū)的土壤有機碳分布。

盡管各種方法所得空間插值圖有一定的相似之處，但各種方法所得成果圖的細部具有較大差別。KYJZ、KDL所得的插值圖最為相似，圖斑較大，平滑效應明顯，表明這兩種方法所得結果基本一致。而KDLTR和KDMTR法所得插值圖明顯較上述兩種方法所得圖斑來得破碎，表明對于研究區(qū)而言，土壤類型信息較土地利用類型信息和地貌類型信息對反映土壤有機碳含量的分布影響較大。KDLTR法所得圖斑最為破碎，分辨率最高，表明其更好地表達了研究區(qū)土壤有機碳含量的空間變異特點，采用該方法具有較好的空間插值效果。

2.3 結合不同類型信息的空間插值模型預測精度分析

ArcGIS軟件所進行的交叉驗證是利用參加空間插值的數據組進行的自我驗證，僅能反映基于該數據組的最佳插值模型及其參數，從理論上說具有一定的“欺騙性”。在ArcGIS軟件所進行的交叉驗證的基礎上，為檢驗不同樣點密度下，各種不同空間插值方法的預測精度，本研究對驗證樣點的預測值與相同位置的實測值進行相關性分析，并作散點圖（圖2）進行檢驗，擬合出兩者間的回歸方程。同時，采用均方根誤差（RMSE）（圖3）來衡量各種方法預測結果的準確性和預測方法的系統誤差。

從圖2中可以看出，各方法的實測值與預測值均達到了極顯著相關水平（p＜0.01），且其擬合線的斜率不同。同時各種格網密度下采用的各種空間插值方法間的均方根誤差最大值與最小值相差均達到0.5 g／kg以上，表明在同一格網密度內采用的各種預測方法所得的預測結果的精度具有較明顯的差距。

在0.5 km×0.5 km的樣點密度中，KDLTR法的預測值和實測值的相關系數最大（r＝0.786＊＊）其擬合線的斜率也較接近于參考線，KDL法中的相關系數最小（r＝0.744＊＊），其擬合線的斜率與參考線的斜率相差也較大，而其它方法的相關系數介于兩者之間，由大到小分別為 KTR、KDMTR、KYJZ和KDM。同時，KDLTR法所得預測結果的均方根誤差最小，為3.641 9 g／kg，KDL法的均方根誤差最大，達4.422 6 g／kg，兩者相差0.780 7 g／kg；其余方法中均方根誤差由小到大依次為KTR、KDMTR、KDM和KYJZ。KDLTR法相對于KYJZ法，其均方根誤差下降了17.38%，即其精度提高了17.38%。從各種方法的均方根誤差還可看出（圖3），KDLTR、KDMTR和KTR這三種結合了土壤類型信息的方法均較其它方法的預測精度有較為明顯的提高，表明在縣級尺度下，土壤類型信息對于土壤有機碳含量的影響較大。

圖2 不同方法的實測值與預測值線性回歸分析

圖3 不同格網密度內各種空間插值方法的均方根誤差

在1 km×1 km的樣點密度中KDMTR法的預測值和實測值的相關系數最大（r＝0.792＊＊），KDM法兩者的相關系數最?。╮＝0.689＊＊）。其它方法的相關系數介于兩者之間，由大到小分別為KTR、KDLTR、KDL和KYJZ。同時，KDLTR法所得預測結果的均方根誤差最小，為4.650 5 g／kg，KDM法的均方根誤差最大，為5.232 0 g／kg，兩者相差0.581 5 g／kg；其余方法中均方根誤差由小到大依次為KTR、KDMTR、KYJZ和KDL。KDLTR法相對于KYJZ法，其均方根誤差下降了9.64%。

在2 km×2 km的樣點密度中KTR法的預測值和實測值的相關系數最大（r＝0.752＊＊），KDM 法兩者的相關系數最?。╮＝0.622＊＊）。其它方法的相關系數介于兩者之間，由大到小分別為KDMTR、KDLTR、KYJZ和KDL。同時，KTR法所得預測結果的均方根誤差最小，為4.680 5 g／kg，KDM法的均方根誤差最大，為5.578 4 g／kg，兩者相差0.897 9 g／kg；其余方法中均方根誤差由小到大依次為KDLTR、KDMTR、KYJZ和KDL。KTR法相對于KYJZ法，其均方根誤差下降了10.04%。

在4 km×4 km的格網密度中，相關系數最大的為KTR法，（r＝0.617＊＊），KDM 法的相關系數最小，為0.461＊＊。同時，KYJZ法所預測結果的均方根誤差最小，為6.101 0 g／kg，各種方法所得均方根誤差均大于6.0 g／kg，平均達6.39 g／kg，比2 km×2 km樣點密度的平均均方根誤差增加了25.72%，其主要原因仍應歸于樣點數量的減少，不能充分表征土壤有機碳的空間變異。研究結果表明4 km×4 km的樣點密度已不適用于南方丘陵區(qū)縣級尺度的土壤有機碳含量的空間預測。

2.4 樣點密度變化對各種預測方法的精度影響

將各種預測方法在不同樣點密度下所得預測結果與實測值的相關系數進行比較分析可知，除KTR和KDMTR法外，其余方法的相關系數均隨格網密度減小而減小。

比較各種預測方法在不同樣點密度下所得預測結果與實測值的均方根誤差（圖4），可知各種預測方法的均方根誤差隨著格網密度的減少而增大。0.5 km×0.5 km的樣點密度下，其各種方法的空間插值精度均較其它密度等級有明顯提高；而1 km×1 km和2 km×2 km兩種樣點密度的空間插值精度差距很小，表明1 km×1 km的樣點密度對其耕層有機碳含量空間預測精度的影響不明顯；4 km×4 km的樣點密度下，各種空間插值方法所得結果的誤差均有較大幅度的提高。就各種方法而言，KTR法在各個樣點密度等級下均有較良好的表現，是一種較好的空間插值方法，也表明在縣級尺度上，土壤類型是影響土壤有機碳的重要因素，而僅結合土地利用現狀或地貌類型信息的克里格插值方法，對于土壤有機碳含量的空間預測效果則不理想。KDLTR法在樣點密度大于1 km×1 km時，其預測結果的精度最高；當樣點密度小于等于1 km×1 km時，則KTR法的預測精度最高。

圖4 不同空間插值方法的精度對格網密度的響應

以上分析表明，當樣點密度小于等于4 km×4 km的格網密度時，由于樣點數過少，不適用于南方丘陵區(qū)縣級尺度的土壤有機碳空間預測。如受條件制約，僅能布設到該密度等級時，采用結合類型信息的克里格插值法的預測精度不如采用原始有機碳含量數據直接進行普通克里格插值法所得結果的預測精度。因此，在縣級尺度上，當研究所需較高的空間預測精度時，樣點的布設須達到0.5 km×0.5 km及以上的樣點密度，此時應采用結合土地利用現狀和土壤類型信息的普通克里格法進行空間預測；當研究僅需達到中等精度的空間預測結果時，則僅需以2 km×2 km的格網密度布設樣點，此時宜采用結合土壤類型信息的普通克里格法進行空間預測。

3 結論

（1）龍海市土壤有機碳含量在縣級研究區(qū)域上存在較強的空間變異性。

（2）結合不同類型信息空間插值模型對于土壤有機碳區(qū)域分布的預測精度均有不同程度的提高，同時樣點密度直接影響了空間插值方法的選擇。根據研究目的的不同，當研究只要求獲得中等精度的空間預測結果時，應按2 km×2 km的格網密度布設調查樣點，并采用結合劃分到土屬級別的土壤類型信息進行普通克里格法（KTR）插值；當研究要求獲得高精度的空間預測結果時，則需按0.5 km×0.5 km及以上的格網密度進行樣點布設，并采用結合土地利用現狀類型和土壤類型信息的普通克里格法（KDLTR）進行空間插值；而格網密度小于4 km×4 km時，由于樣點數太少，不能得出準確反映土壤有機碳含量的空間分布，若受條件限制，僅能得到該格網密度時，應直接采用普通克里格法（KYJZ）進行空間預測。

（3）在南方丘陵區(qū)縣級尺度上，土壤類型較土地利用類型和地貌類型對土壤有機碳含量的空間分布影響大。凡是涉及結合土壤類型信息的克里格方法（KTR、KDLTR、KDMTR）所得空間插值結果均明顯好于未結合土壤類型信息的方法（KYJZ、KDL、KDM）。因此在縣級尺度上研究土壤有機碳的空間分布時最好考慮土壤類型的影響。

［1］華孟，王堅.土壤物理學［M］.北京：北京農業(yè)大學出版社，1992：214-243.

［2］楊玉玲，文啟凱，田長彥.土壤空間變異研究現狀及展望［J］.干旱地區(qū)農業(yè)研究，2001，18（2）：50-55.

［3］胡克林，李保國，林啟美，等.農田土壤養(yǎng)分的空間變異性特征［J］.農業(yè)工程學報，1999，15（3）：33-38.

［4］王軍，邱揚.土地質量的空間變異與尺度效應研究進展［J］.地理科學進展，2005，24（4）：28-35.

［5］張淑娟，方慧，何勇.精細農業(yè)田間信息采樣策略［J］.農業(yè)機械學報，2004，35（4）：88-92.

［6］Thomas J S，David W M.The Spatial variation of plantavailable phosphorus in pastures with contrasting management［J］.Soil Science Society of America Journal，2003，67：826-836.

［7］Wu J，Norvell W A，Hopkins D G，et al.Improved prediction and mapping of soil copper by kriging with auxiliary data for cation-exchange capacity ［J］.Soil Science Society of America journal，2003，67（3）：919-927.

［8］姜城，楊俐蘋.土壤養(yǎng)分變異與合理取樣數量［J］.植物營養(yǎng)與肥料學報，2001，7（3）：262-270.

［9］王宏庭，金繼運，王斌，等.土壤速效養(yǎng)分空間變異研究［J］.植物營養(yǎng)與肥料學報，2004，10（4）：249-354.

［10］Wang X J，Qi R.The effect of sampling design on spatial structure analysis of contaminated soil［J］.The science of the total environment，1998，224：29-41.

［11］Liu T L，Juang K W，Lee D Y.Interpolating soil properties using kriging combined with categorical information of soil maps［J］.Soil Science Society of America Journal，2006，70（4）：1200-1209.