基于GIS與分區Kriging的采煤沉陷區土壤有機碳含量空間預測

徐占軍，張 媛，張紹良，李樂樂，余明成

（1. 山西農業大學資源環境學院；晉中 030801；2. 中國礦業大學環境與測繪學院；徐州 221008）

0 引 言1

2014年，中國煤炭產量為2.6×109t，占世界煤炭總產量的 46.9%[1]。煤炭開采過程中的一系列環節占用甚至損毀了大量的土地資源，導致地上植被被移除，土壤動物和其他動物因生境喪失而消亡，礦區的農作物減產，生態環境遭到嚴重破壞[2]。而且，煤礦開采塌陷的土地大部分是平原農業區的可耕地，即華北平原、東北平原和長江中下游平原的農業區[3]。相關研究表明，中國年塌陷的耕地面積約為200 km2[4]。人類煤炭開采活動對農業生態系統危害嚴重[5]，對農田土壤有機碳碳庫擾動十分劇烈[6]。

由于農田的土壤有機碳庫是減少陸地生態系統碳排放的最大潛在因素[7-8]，中國以及世界上的其他煤炭開采大國必須在區域尺度上就煤炭開采對區域農田碳儲量的影響進行定量的研究，才可以更好地對煤炭開采區的土壤有機碳庫進行科學管理，同時也實現煤炭低碳開采。而適宜于煤炭開采沉陷區的區域土壤有機碳含量空間預測模型是必須要解決的科學問題。但是由于礦區受到人類影響較大，區域內部的土壤都不同程度地出現了地表破壞、土壤侵蝕、地面沉陷、植被毀損等問題[9-13]，無論哪一方面的變化都會對土壤碳庫產生影響，使礦區土壤有機碳庫通常存在比較強的空間變異性[14-19]，因此，有必要結合煤炭開采活動引起的擾動影響因素對礦區土壤有機碳庫進行估算，以更好地了解煤炭開采活動帶來的區域尺度上土壤有機碳庫含量的變化。

當前關于區域土壤含量預測主要有 2種思路，一種是基于遙感數據對區域土壤有機碳庫預測，另一種是運用地統計學方法進行預測。但目前，基于遙感的方法存在著一個主要的缺陷，即尚未建立航天遙感與土壤有機碳庫之間的定量關系，因此地統計學思路就成為預測區域土壤有機碳庫的主要方法。

Kriging空間插值法是應用最廣泛的土壤屬性空間預測方法[20-23]。但是普通Kriging比較適用于地質及人文環境等各方面條件相對比較均質的區域，在空間區域較小但地形復雜，土壤屬性受人類活動擾動、變化特別強烈的地區，其應用效果并不理想。

隨著計算機技術的不斷進步，結合輔助信息的Kriging方法已經得到了廣泛的應用，該方法可以提高預測精度。已有研究表明，許多變量可以作為Kriging插值法的輔助變量，提高土壤有機碳庫空間預測的精度[24]。例如Chai等[25]以北京平谷區為研究區，分別利用基于約束性極大似然模型的最佳線性無偏估計法和回歸Kriging法，對區域土壤有機質含量進行了空間預測，并比較了二者的預測精度；Mishra 等[20]以印第安納州為研究區，分別使用剖面深度分布函數和普通Kriging模型，預測了表層土壤有機碳的含量，并比較了 2種預測模型的預測精度。Aladamat R.等[26]建立了全球土壤有機碳建模系統，并結合研究區土地利用的變化情況，預測了約旦全國尺度上2000年到2030年的土壤有機碳的儲量。

本文提出在煤炭開采劇烈擾動的開采沉陷區，利用結合輔助變量的 Kriging法來預測土壤有機碳含量的思路。在高潛水位煤礦區，煤炭開采沉陷形成了不同的積水區，根據不同的積水情況人們采取了不同的土地利用措施，導致不同區域之間土壤有機碳含量呈現出很大的差異。本文以徐州九里煤炭開采沉陷區為研究區，采用分區Kriging法，以沉陷積水情況為分區輔助變量，進行了土壤有機碳含量預測；同時也利用普通Kriging法對沉陷區土壤有機碳庫含量進行了預測。通過對 2種方法的預測精度的比較，選出精度更好的礦區土壤有機碳庫儲量估算模型，作為礦區煤炭開采沉陷區內土壤有機碳庫含量空間預測的方法。

1 材料與方法

1.1 研究區概況

九里采煤塌陷區（34°13'39''–34°26'16''N, 117°06'21''–117°12'16''E）位于江蘇省徐州市九里區，面積約為42.15 km2，是徐州市張小樓煤礦、龐莊煤礦、夾河煤礦的煤炭開采沉陷區。3礦區的井田邊界以及沉陷區域如圖1所示（紅線范圍內）。研究區煤礦的開采方式為井工開采，選取長壁采煤的技術。這種采煤技術造成了地面塌陷，并伴隨著一系列的生態環境問題，主要包括：（1）煤炭開采導致農田大面積塌陷；（2）煤炭開采活動導致一些位于高潛水位的農田變成了沉陷積水水域。

圖1 研究區位置以及采樣點分布圖Fig.1 Location of study area and distribution of soil sampling sites

研究區屬于中國華東平原區，平均海拔在 35～44 m之間，地下水位較淺，屬于溫暖季風氣候區，年平均氣溫14 ℃，年平均無霜期為200～220 d，年平均降水量為800～-930 mm。該區是古黃河沖積平原，水稻土、砂姜土為該區主要土壤類型。農業生產按照冬小麥/夏玉米的方式進行輪作，屬于中國東部的糧食生產區。

1.2 土壤樣品采集與處理

采樣點包括預測樣點和驗證樣點 2部分，預測樣點采樣采用的方法是正方形格網法。根據研究區的實際面積，格網大小為1000 m 1000 m×，共采集了54個預測樣點。本研究只采用了每個土樣的0～20 cm土壤有機碳含量數據。在采樣過程中，用手持GPS機記錄下每一個采樣點的地理坐標，其中沉陷濕地區內采集8個土壤樣點，季節性積水區內采集13個土壤樣點，沉陷未積水區內采集33個土壤樣點。另外，在整個沉陷區內再隨機采集18個樣點作為驗證樣點，與 2種估算模型的預測值對比，對模型所得出的結果進行精度評價。外業采樣完成后對土壤有機碳含量采用重鉻酸鉀氧化法進行測定[27]。本次采樣工作是在2010年10月份研究區域內所有的農作物收割完成以后進行的。

1.3 煤炭開采沉陷區土壤有機碳含量空間預測方法

對區域土壤理化性質進行空間預測常用的方法是Kriging空間插值法。它是基于變異函數模型在有限區域內對區域化變量的取值進行無偏最優估計的空間插值方法。其基本原理和方法在許多文獻中均有詳細描述[23,28-29]，是基于區域化變量理論，利用原始數據和半方差函數的結構性，計算可獲取變量的線性加權組合，從而對待估點預測值進行無偏的最佳估計，而且這種方法能夠計算估計值的可靠性程度（估計值方差）。

研究區包括多個煤礦，且處于高潛水位區。由于每個礦區對煤炭資源的開采程度不同，導致各個礦區內部的地表沉陷程度不同，進而形成了不同程度的沉陷積水區（包括沉陷未積水區、季節性積水區和常年積水濕地區），人們針對礦區不同的積水情況采取了不同的土地利用方式，以及農業耕作措施，因而土壤有機碳含量受到影響也會產生一定的差異[30]。在高潛水位礦區由開采沉陷導致的沉陷積水情況，是影響沉陷區土壤有機碳庫的主導因素，因此根據不同的積水情況來分區進行Kriging插值，將樣本中土壤有機碳庫的含量分為 2部分，即各個積水區的均值和殘差。其計算公式如下：

式中 Z （ Xkj）表示各個采樣點的有機碳含量， Xkj表示的是土樣所在的位置， u（ tk）表示各個積水區有機碳含量的均值， tk表示的是土樣所在沉陷區的積水情況：主要分為常年積水濕地區、季節性積水區以及未積水區 3個部分。其中，k表示不同的區域積水情況，j表示不同的土壤樣點。r(xkj)表示殘差，即各個采樣點的有機碳減去有機碳均值的差。其中，均值間的差異反映的是沉陷積水情況不同區域中有機碳含量的變異性，而殘差之間的差異反映的是沉陷積水情況相同的地區有機碳含量的變異性。在這種方法中，將殘差看作是新的區域變量來進行Kriging插值，其變異函數 γr（ h）及待估點Xkj的預測公式分別為式（2）和式（3）。每一個待估點的土壤有機碳庫含量預測值 Z*（xkj）為所屬沉陷積水情況區的均值u(tk)與

所估計的殘差r*（xkj）之和（式4）。

式中h為2個樣本間的點空間分隔距離，N(h)是分隔距離為h時的樣本點對總數。

式中kjλ為權重系數。

在本次研究中，首先用SPSS來對土壤有機碳庫數據及其殘差進行統計分析，然后利用ArcGIS 10.0軟件中的Geostatistical Analyst 模塊，檢驗數據的質量和分布，本次的采樣數據滿足正態分布和隨機性檢驗，用該軟件進行空間插值并生成空間插值圖。

1.4 預測結果準確性檢驗

本次對于研究結果準確性的檢驗是通過比較18個驗證樣點位置處的預測值與采樣實測值之間的相關系數（r）[31]，以及Isaaks和Shrivastava定義的3個驗證指標，即平均誤差（mean error，ME）、平均絕對誤差（mean absolute error，MAE）和均方根誤差（root mean square error，RMSE）來評估不同插值方法的性能[32]。其中，相關系數（r）這一指標是用以反映變量之間相關關系密切程度的。均方根誤差RMSE實際上是觀測值與真值偏差的平方和再與觀測次數N來求比值，然后再求平方根來計算得到的。

式（5）～（7）中，N表示驗證點的總數量，本次研究中N為18，xoi表示驗證點的實測值，xpi表示根據估算模型得到的驗證點預測值。

在實際的測量過程中，由于對樣點進行的觀測次數總是有限的，因此只能用最可信賴的值來代替真值，而對一組測量中的特大或特小誤差，標準誤差的反映非常敏感，所以，可以用標準誤差來很好地反映出預測的精確度。r值越大、RMSE值越小，說明所預測精度越高, 反之則說明預測精度越低。

2 結果與分析

2.1 土壤有機碳的統計分析

首先對試驗區所采集的土壤樣本進行了描述性統計分析（表 1）。分析顯示，區域內采樣點的土壤有機碳含量符合正態分布規律，不同積水區內土壤有機碳庫含量的差異較為顯著，因此能夠進行Kriging插值；同時也為下一步對土壤有機碳庫含量進行半方差函數結構分析以及相應的分區Kriging插值方案的確定提供了依據，并且減少了統計分析過程中所產生的不確定性。根據表 2可知，3個區域內土壤有機碳含量從小到大依次為未積水區、季節性積水區、濕地，這說明土壤有機碳庫的含量在不同的積水區域內有很大的不同。而且積水區內土壤有機碳庫的變異系數為 0.07，比未積水區內土壤有機碳庫的變異系數0.13要小，尤其在季節性積水區內其標準差為 0.86，也較小，說明積水區內的土壤有機碳庫含量離散程度較小，比較集中，這也從另一方面驗證了將沉陷區內積水量的變化作為輔助變量的分區Kriging插值法的可行性。

表1 沉陷區域外業采樣樣本土壤有機碳庫含量描述性統計結果Table 1 Descriptive statistical results of soil organic carbon content of samples in subsidence area

2.2 地統計分析及Kriging插值

利用ArcGIS 10.0，對采樣點樣本土壤有機碳含量的原始數據、去除區域積水分區內有機碳含量均值以后的殘差半方差函數及其擬合參數進行分析[33-36]，其結果見表3。

表2 沉陷區不同積水情況下土壤樣本有機碳描述性統計結果Table 2 Descriptive statistics results of soil organic carbon samples with different water accumulation in subsidence area

表3 土壤有機碳含量及其殘差半方差函數模型和參數Table 3 Soil organic carbon content and its residual semi variance function model and parameters

表 3顯示了樣本有機碳含量原始數據和土壤有機碳殘差各自的半方差函數擬合指數和球狀模型的參數。根據塊金值的相關定義，理論上當采樣點的距離為 0時，半變異函數值應為 0，但由于存在測量誤差以及空間變異，使得 2個采樣點非常接近時，它們的半變異函數值不為0，即存在塊金值，因此塊金值反映的是最小抽樣尺度以下變量的變異性及測量誤差。結合對表 3的分析，可以看出采樣點土壤有機碳含量與殘差的塊金值分別為0.91和 0.01，土壤有機碳的空間異質性較大，殘差的空間異質性比較小，說明該地區不同積水沉陷區域內部由于試驗誤差和一些隨機因素引起的空間變異都很小。

C是結構方差，表示系統屬性或區域變量最大空間變異，是由土壤母質、地形、氣候等結構性因素引起的變異。由于該沉陷區域開采前后氣候并未發生變化，土壤有機碳庫的空間變異性基本是由于煤炭開采造成的地表形變和沉陷積水等結構因素引起的變異，該沉陷區域土壤有機碳結構方差C是15.17，這種結構性因素是該區域土壤有機碳庫變異的主導因素。

C/(C0+C)即塊金值與基臺值的比值，被稱為空間相關度，表示可度量空間自相關的變異所占的比例，該值的大小表示系統變量的空間相關性的程度。由表3可知2種方法得出的 C/(C0+C)分別為 0.998與 0.943，均大于75%，這說明其相關性很強，尤其是樣本的土壤有機碳含量的相關性極強，達到了0.998。

根據表 3可知，與原始數據相比，土壤有機碳殘差擬合函數的塊金值很小，基臺值降低，變程也降低，這主要是因為，殘差是土壤有機碳含量減去各個積水區內的有機碳含量均值得到的，而不同積水區域內部的土壤有機碳含量的變化比較小，相對比較均勻，這也證明了對沉陷區內土壤有機碳含量進行預測時考慮不同積水沉陷區的合理性，間接地說明了以分區對沉陷區土壤碳庫儲量進行Kriging插值預測的有效性和可行性。

2.3 兩種方法預測精度的對比分析

采用研究區18個驗證點處的土壤有機碳含量實測值與模型的預測值進行對比，對分區Kriging模型的預測精度進行評價。圖2是預測值與實測值的散點分布圖，由2種方法所得到的線性回歸方程顯示，以區域積水情況為輔助變量進行的Kriging插值得到的預測值與實測值的決定系數（0.756 4）明顯高于普通的Kriging插值法（0.508 6）。

表4為2種預測方法的驗證點處土壤有機碳含量預測的準確性評估指標。ME表示預測的平均偏差，從表4中可以看出，普通Kriging（0.020 2 g/kg）和結合積水情況的Kriging（0.021 8 g/kg）的ME均接近0，說明2種方法預測的平均偏差都比較小，總體預測精度都比較高。平均絕對誤差（MEA）反映了預測值誤差的實際情況，由表4可知，普通Kriging的MAE為1.851 1 g/kg，結合積水情況的Kriging的MAE為1.287 8 g/kg，因此，結合積水情況的Kriging預測實際誤差更小。作為平均值和方差偏差的聯合度量指標，如表4所示，直接進行Kriging插值對土壤有機碳含量的預測RMSE為0.55，明顯高于分區Kriging插值（0.35）。綜上所述，考慮區域內土壤積水情況的Kriging插值法對土壤有機碳含量的預測值精度更高。

圖2 兩種方法得到的土壤有機碳預測值與實測值散點圖Fig.2 Scatter diagram of predicted and measured values of soil organic carbon by 2 methods

表4 兩種方法預測土壤有機碳含量的準確性評估指標Table 4 Accuracy assessment indices of estimated soil organic carbon content with 2 methods

2.4 土壤有機碳的空間分布特征

通過 2種方法預測得到的土壤碳庫儲量空間分布如圖3所示，可以發現沉陷區2種土壤有機碳含量空間預測方法（普通Kriging和分區Kriging）得到的土壤有機碳殘差的空間分布格局，而分區Kriging預測的土壤有機碳含量圖的空間分布特征及空間遞變規律更加明顯，普通Kriging得到的土壤有機碳含量范圍為：9.34～16.252 g/kg，而分區Kriging插值的結果為9.338～18.058 g/kg，2種方法的估計范圍大致相同，總體上看研究區土壤有機碳分布的空間格局基本一致，有機碳含量由中部向四周逐漸減少，中部的土壤有機碳含量最高，據實地調查該區域屬于濕地區，與其他土地利用相比，濕地土壤長期飽和度較低，其凋落物和有機質分解速度較慢，這增加了濕地土壤中的土壤有機碳含量[37]。

圖3 兩種方法得到的土壤有機碳含量空間分布Fig.3 Spatial distribution of soil organic carbon content by 2 methods

西南部的土壤有機碳含量較少，該區域屬于未積水區，由于礦區煤炭開采對農田土壤和植被擾動影響嚴重，造成植被NPP (net primary productivity) 降低，從而減少了土壤有機碳的一個重要補給來源，而且開采沉陷形成很多沉陷坡面，土壤有機碳由于土壤侵蝕的原因流失嚴重[38-39]。從整體上來說 2種空間插值方法都能夠反映出整個研究區的土壤有機碳分布情況。

但是在局部分布上，2種預測方法之間存在著一定的差異。對比分析發現，直接進行Kriging插值得到的土壤有機碳空間分布的圖斑比較連續規整，這明顯跟研究區的實際情況不太相符，說明平滑效應使不同積水情況區域之間的差異性降低，其差值結果只能反映研究區土壤有機碳的大致分布情況，卻不能更為精確地反映研究區的實際情況，比如圖3a局部S1區以及S2區域由于個別預測樣點的土壤有機碳含量偏高，造成區域內土壤有機碳含量空間插值結果異常偏高。普通Kriging插值是基于采樣點的空間分布特征以及樣點與樣點之間的空間位置關系來進行空間插值，事先未對插值樣點進行分類，所以無法消除鄰近的不同類別空間預測樣點對空間插值結果的影響。而分區Kriging由于按照影響沉陷區土壤有機碳含量的最大因素——沉陷積水情況進行分區，就可以有效消除這方面的誤差影響。而且，與分區Kriging的模擬結果相比，直接Kriging法模擬的礦區土壤有機碳分布特征規律性不明顯，忽略了土壤有機碳空間平穩過渡。因此，分區Kriging模擬更能反應土壤有機碳的空間遞變特點，有利于分析不同因素對土壤有機碳空間分布的影響，能更好地反映研究區域內部土壤有機碳含量的分布情況。

分區Kriging空間插值的思想在相關領域也得到了學者的認可，Zhang等[31]通過對比在土壤現場采樣過程中分配樣點的 4種不同模式，即不分類的網格模式，基于土壤類型的分配模式，基于土地利用類型的分配模式以及基于土地利用模式-土壤類型的分配模式下得到的土壤有機碳的變異系數，得出評估土壤有機碳空間分布的最有效方法是基于土地利用模式-土壤類型分配模式。同樣地，Sandeep等[40]在研究美國賓夕法尼亞州的土壤有機碳時，采用結合環境變量的地理加權回歸克里金法，與原來的回歸克里金法進行對比發現，因為前者考慮了空間非平穩性以及殘差的空間自相關性，其精度會提高。所以在受人類煤炭開采活動擾動劇烈的沉陷區土壤有機碳含量的預測應用中，其預測精度也高于回歸克里金法。

3 結 論

本文提出了利用分區 Kriging法來對煤炭開采沉陷區土壤有機碳含量進行空間插值預測的方法，并將該方法與傳統的直接Kriging方法得到的預測結果進行對比，發現：

1）結合區域內部積水情況來進行的分區 Kriging法得到的研究區土壤有機碳含量的范圍為：9.338～18.058 g/kg，而直接Kriging方法得到的范圍為：9.34～16.252 g/kg。

2）從區域預測精度上分析，與直接進行 Kriging插值相比，分區Kriging預測方法精度更高。研究區內由于煤炭開采造成了不同的沉陷積水區，土地利用方式也各不相同，導致區域內土壤有機碳含量出現差異。結合區域積水情況的分區Kriging，可以消除鄰近不同積水情況區的采樣點土壤有機碳含量對其空間插值預測精度的影響，提高土壤有機碳含量的預測精度。

綜上所述，可以選擇結合區域積水情況的分區Kriging空間插值作為煤炭開采沉陷區土壤有機碳含量空間預測的預測方法。在煤礦區，優化的土壤有機碳空間預測方法，也有利于更精確地進行土壤有機碳動態模擬和理解本地土壤有機碳庫的時空演化，為礦區土地低碳復墾乃至區域內的土地資源的低碳利用提供更加科學的依據。

[1] BP. The BP Statistical Review of World Energy 2014[EB/OL]. http://www.bp.com/statisticalreview, 2018-03-12.

[2] 原野, 趙中秋, 白中科, 等. 露天煤礦復墾生物多樣性恢復技術體系與方法：以平朔礦排土場為例[J]. 中國礦業,2017, 26(8): 93－98.Yuan Ye, Zhao Zhongqiu, Bai Zhongke, et al. Technology system and method of biodiversity restoration for the reclamation of opencast coal mine:a case study from the dumps in Pingshuo mine[J]. China Mining Magazine, 2017,26(8): 93－98. (in Chinese with English abstract)

[3] 羅愛武. 淮北市采煤塌陷區土地復墾研究[J]. 安徽師范大學學報（自然科學版）, 2002, 25(3): 286－289.Luo Aiwu. The land restoration in the sunk areas of coal extraction in HuaiBei City[J]. Journal of AnHui Normal University (Nature Science), 2002, 25(3): 286－289. (in Chinese with English abstract)

[4] 彭蘇萍, 趙建慶. 中國西部煤礦區生態環境控制及改善研究[C]//中國科協2000年學術年會. 2000.

[5] 侯湖平, 徐占軍, 張紹良, 等. 煤炭開采對區域農田植被碳庫儲量的影響評價[J]. 農業工程學報, 2014, 30(5): 1－9.Hou Huping, Xu Zhanjun, Zhang Shaoliang, et al. Effect evaluation on vegetation carbon pool of region agro-ecosystem by coal mining in mining area[J].Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(5): 1－9.(in Chinese with English abstract)

[6] 王坤. 高潛水位采煤塌陷區充填復墾土壤碳動態研究[D].徐州：中國礦業大學, 2016.Wang kun. Study on Carbon Dynamics of Reclaimed Soil in Coal Mining Subsidence Area with High Groundwater Level[D]. Xuzhou: China University of Mining and Technology, 2016.

[7] Miller B A, Koszinski S, Wehrhan M, et al. Comparison of spatial association approaches for landscape mapping of soil organic carbon stocks[J]. Soil, 2015(1): 217－233.

[8] Zhang G S, Ni Z W. Winter tillage impacts on soil organic carbon, aggregation and CO2emission in a rainfed vegetable cropping system of the mid–Yunnan plateau, China[J]. Soil& Tillage Research, 2017, 165: 294－301.

[9] Lan C Y, Shu W S, Wong M H. Revegetation of lead/zinc mine tailings at shaoguan, guandong province, China:Phytotoxicity of the tailings[J]. Studies in Environmental Science, 1997, 66(97): 119－130.

[10] Yan Chu. The influence of coal mining on environmental quality of mining area and the research progress of repairing technology[J]. Advances in Environmental Protection, 2016,6(1): 1－6.

[11] Anthony B. Restoration of mined land-using natural processes[J]. Ecological Engineering, 1997, 8(4): 255－269.

[12] Zhang Ling, Wang Jinman, Bai Zhongke, et al. Effects of vegetation on runoff and soil erosion on reclaimed land in an opencast coal-mine dump in a loess area[J]. Catena, 2015,128(5): 44－53.

[13] Huang Yi, Tian Feng, Wang Yunjia, et al. Effect of coal mining on vegetation disturbance and associated carbon loss[J]. Environmental Earth Sciences, 2013, 73(5): 1－14.

[14] álvaro-Fuentes J, Easter M, Paustian K. Climate change effects on organic carbon storage in agricultural soils of northeastern Spain[J]. Agriculture Ecosystems &Environment, 2012, 155(155): 87－94.

[15] Yang X M, Drury C F, Wander M M, et al. Evaluating the effect of tillage on carbon sequestration using the minimum detectable difference concept[J]. Pedosphere, 2008, 18(4):421－430.

[16] Mishra U, Ussiri David A N, Lal R. Tillage effects on soil organic carbon storage and dynamics in Corn Belt of Ohio USA[J]. Soil & Tillage Research, 2010, 107(2): 88－96.

[17] Rosemary F, Vitharana U W A, Indraratne S P, et al.Exploring the spatial variability of soil properties in an Alfisol soil catena[J]. Catena, 2017, 150(3): 53－61.

[18] Teng Mingjun, Zeng Lixiong, Xiao Wenfa, et al. Spatial variability of soil organic carbon in Three Gorges Reservoir area, China[J]. Science of the Total Environment, 2017,599(6): 1308－1316.

[19] Wu Lizhi, Li Long, Yao Yunfeng, et al. Spatial distribution of soil organic carbon and its influencing factors at different soil depths in a semiarid region of China[J]. Environmental Earth Sciences, 2017, 76 (19): 654.

[20] Mishra U, Lal R, Slater B, et al. Predicting soil organic carbon stock using profile depth distribution functions and ordinary kriging[J]. Soil Science Society of America Journal,2009, 73(2): 614－621.

[21] Liu Yaolin, Guo Long, Jiang Qinghu, et al. Comparing geospatial techniques to predict SOC stocks[J]. Soil &Tillage Research, 2015, 148(3): 46－58.

[22] 杜挺, 楊聯安, 張泉, 等. 縣域土壤養分協同克里格和普通克里格空間插值預測比較——以陜西省藍田縣為例[J].陜西師范大學學報(自然科學版), 2013, 41(4): 85－89.Du Ting, Yang Lian'an, Zhang Quan, et al. Spatial prediction comparison of soil nutrient between ordinary Kriging and Cokring at county scale-A case study in Lantian county of Shaanxi province[J]. Journal of Shaanxi Normal University(Natural Science Edition), 2013, 41(4): 85－89. (in Chinese with English abstract)

[23] 施周, 閆杭召, 畢晨, 等. 基于地統計學——克里格插值法的村鎮地表水體水質監測[J]. 環境工程學報, 2017(4):2607－2613.Shi Zhou, Yan Hangzhao, Bi Chen, et al. Surface water quality monitoring in town based on Geotatistics-Kriging interpolation[J]. Chinese Journal of Environmental Engineering, 2017(4): 2607－2613. (in Chinese with English abstract)

[24] Dai Fuqiang, Zhou Qigang, Lv Zhiqiang, et al. Spatial prediction of soil organic matter content integrating artificial neural network and ordinary kriging in Tibetan Plateau[J].Ecological Indicators, 2014, 45(5): 184－194.

[25] Chai Xurong, Shen Chongyang, Yuan Xiaoyong, et al.Spatial prediction of soil organic matter in the presence of different external trends with REML-EBLUP[J]. Geoderma,2008, 148(2): 159－166.

[26] Aladamat R, Rawajfih Z, Easter M, et al. Predicted soil organic carbon stocks and changes in Jordan between 2000 and 2030 made using the GEF SOC Modelling System[J]. Agriculture Ecosystems & Environment, 2007, 122(1): 35－45.

[27] Fang Yu, Yan Zhilei, Chen Jichen, et al. Effect of chemical fertilization and green manure on the abundance and community structure of ammonia oxidizers in a paddy soil[J]. Chilean Journal of Agricultural Research, 2015, 75(4): 487－496.

[28] Wang L X, Okin G S, Caylor K K, et al. Spatial heterogeneity and sources of soil carbon in southern African savannas[J]. Geoderma, 2009, 149(3-4): 402－408.

[29] Sun Zhili, Wang Jian, Li Rui, et al. LIF: A new Kriging based learning function and its application to structural reliability analysis[J]. Reliability Engineering & System Safety, 2017, 157(4): 152－165.

[30] Shrestha R K, Lal R. Changes in physical and chemical properties of soil after surface mining and reclamation[J].Geoderma, 2011, 161(3–4): 168－176.

[31] Zhang Z Q, Yu D S, Shi X Z, et al. Effect of sampling classification patterns on SOC variability in the red soil region, China[J]. Soil & Tillage Research, 2010, 110(1): 2－7.

[32] Isaaks E H, Srivastava R M. An Introduction to Applied Geostatistics[M]. Oxford University Press, 1989, 33(33): 483－485.

[33] Lévesque J, King D J. Airborne digital camera image semivariance for evaluation of forest structural damage at an acid mine site[J]. Remote Sensing of Environment, 1999,68(2): 112－124.

[34] Li Xinrong. Influence of variation of soil spatial heterogeneity on vegetation restoration[J]. China Science:Earth Science, 2005, 48(11): 2020－2031.

[35] 郭凌俐, 王金滿, 白中科, 等. 黃土區露天煤礦排土場復墾初期土壤顆粒組成空間變異分析[J]. 中國礦業, 2015,24(2): 52－59.Guo Lingli, Wang Jinman, Bai Zhongke, et al. Analysis of spatial variability of soil granules in early stage of reclamation at opencast coal minedump in loess area[J].China Mining Magazine, 2015, 24(2): 52－59. (in Chinese with English abstract)

[36] 代富強, 周啟剛, 劉剛才. 基于回歸克里格和遙感的紫色土區土壤有機質含量空間預測[J]. 土壤通報, 2014, 45(3):562－567.Dai Fuqiang, Zhou Qigang, Liu Gangcai. Spatial prediction of soil organic matter contents in a purplish soil region with regression Kriging and remote sensing[J]. Chinese Journal of Soil Science, 2014, 45(3): 562－567. (in Chinese with English abstract)

[37] Qualls R G, Richardson C J. Decomposition of Litter and Peat in the Everglades: The Influence of P Concentrations[M].Springer New York, 2008.

[38] Hu Zhenqi, Xiao Wu. Optimization of concurrent mining and reclamation plans for single coal seam: A case study in northern Anhui, China[J]. Environmental Earth Sciences,2013, 68(5): 1247－1254.

[39] Indorante S J, Jansen I J, Boast C W. Surface mining and reclamation: Initial changes in soil character[J]. Journal of Soil & Water Conservation, 1981, 36(6): 347－351.

[40] Sandeep K, Rattan L, Liu Desheng. A geographically weighted regression kriging approach for mapping soil organic carbon stock[J]. Geoderma, 2012, 189–190 (6): 627－634.