基于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.