基于时空变异的旱地土壤有机碳高效采样策略研究*

姚彩燕，刘绍贵，乔 婷，龙 军，于东升，史学正，邢世和，陈瀚阅，张黎明†

（1. 福建农林大学资源与环境学院，福州 350002；2. 福建农林大学公共管理学院，福州 350002；3. 土壤生态系统健康与调控福建省高校重点实验室，福州 350002；4. 扬州市农业环境监测站，江苏 扬州 225603；5. 土壤与农业可持续发展国家重点实验室（中国科学院南京土壤研究所），南京 210008）

土壤有机碳（Soil Organic Carbon，SOC）作为全球陆地生态系统最大的有机碳库，不仅是评价土壤肥力的重要指标，而且其含量的轻微变化也可导致大气二氧化碳（CO2）浓度显著变化，进而影响全球碳平衡[1]。受气候、地形、土地利用和土壤类型等多方面因素综合影响，SOC具有很高的时空变异性，不同时间和空间位置的SOC含量特征复杂且不确定，准确了解SOC时空分布变化，对于制定科学合理的碳汇管理措施具有重要意义[2]。样点是土壤特性时空变异的信息载体，是土壤调查与空间制图的基础。在进行土壤采样时，任何采样过程都不可能获得土壤的完全信息，只能对一定数量的土壤样本数据通过插值形式实现土壤属性空间分布的定量表达[3]。然而，样点数量过小，土壤样品的偶然性增加，缺乏代表性，空间预测将会产生不可接受的误差；但盲目增大采样数量，不仅预测误差不会持续减小，反而加大了采样工作量和样品分析周期，浪费人力财力物力[4]。因此，明确不同样点数量与SOC预测精度的响应关系是制定高效采样策略的基础。

近年来，国内外学者对不同区域的合理采样数目开展了一系列研究[5-7]。但当前最佳土壤采样数目的确定主要集中在若干个样点数量的土壤属性空间预测精度对比上，并没有量化预测精度与采样数量之间的关系。此外，目前的大多数研究一般局限在单一时间尺度下，未考虑不同时期土壤属性差异对适宜采样数目的影响。众所周知，土壤是动态的，经过一定时间，SOC含量通常会出现一定范围的变化，在某些情况下，SOC时间变异将会超过空间变异的影响[8]。因此，为了准确估计SOC含量，基于样点的预测制图必须要考虑空间和时间维度。旱地作为我国主要的农田类型，占全国耕地面积的76.9%，其SOC密度显著低于世界平均水平，固碳潜力巨大[9]。但由于人类活动以及农田管理措施的影响，该土地利用类型下的土壤碳库呈现强烈的时空变异性，明确不同时期旱地SOC预测精度与采样数量之间的关系是合理制定我国农田管理措施的基础。江苏省北部地区85%的耕地为旱地土壤，尽管面积仅占全国旱地土壤总面积的3.7%，却固定了全国10.2%的旱地SOC，具有较大的“碳汇”空间[10]。此外，该地区旱地的主要土壤类型为潮土、盐土和砂姜黑土，分别占全区旱地土壤总面积的52.67%、24.31%和8.13%，这与黄淮海平原主要旱地土壤类型面积比例相一致，有利于研究成果向大区域拓展[11]。基于此，本研究以江苏省北部3.93×104km2旱地（简称“苏北旱地”）为研究区，利用1980年全国第二次土壤普查和2008年农业部测土配方施肥项目的实测样点数据，通过随机抽样的方式设置20组不同数量的样点数据集，分析这两个时期研究区SOC的最佳采样数量，并在此基础上揭示土壤属性时空变异下预测精度与样点数量之间量化关系，研究结果可为苏北旱地和未来我国第三次土壤普查制定合理的有机碳野外调查方案提供理论依据。

1 材料与方法

1.1 研究区概况

苏北旱地（116°21′～120°54′E，32°43′～35°7′ N）地处黄淮平原与江淮平原的过渡地带，是中国沿海经济带的重要组成部分（图1）。该地区属于亚热带向暖温带过渡的季风性气候，年平均降水量800～1 200 mm，年平均气温13～16 ℃，全年平均日照时数2 000～2 600 h，无霜期约220 d。境内地势平坦，呈西北高、东南低的趋势，地形主要以平原为主，海拔高度在0～50 m之间。研究区面积最大的旱地土壤类型为潮土和盐土，其次为褐土、砂姜黑土和棕壤；而石质土、紫色土和石灰土面积较小。成土母质主要包括湖相冲积物、黄泛冲积物、河海相冲积物、其他各河流冲积物以及下蜀黄土等。苏北是我国重要的粮食产区，种植制度主要是夏玉米和冬小麦轮作[12]。

1.2 数据来源与处理

1980年和2008年表层样点数据（0～20 cm）分别来自全国第二次土壤普查中苏北地区29个县（市、区）在《县级土种志》、《地市级土种志》、《省级土种志》和《中国土种志》中记录的旱地剖面资料及农业部测土配方施肥项目。所有土壤样点信息均包括地理位置描述和理化性质等，其中SOC含量可由有机质含量乘以Bemmelen转换系数（0.58）计算。考虑到样点采集与化学分析过程的不确定性，对两期样点的SOC含量进行正态分布检验后进行异常值剔除处理，所采用的剔除方式是结合数据的常规统计学特征和空间位置，将每个样点的SOC含量与总体及其邻近8个样点的均值和标准差进行比较，如果样点SOC含量在总体均值的五倍标准差之外，且大于或小于邻近样点均值的三倍标准差，则视为异常值剔除，否则予以保留[13]。最后，1980年和2008年可利用的数据集分别为884个和1 081个（图2）。

鉴于矢量数据库可以准确表征SOC预测结果的空间分布特征，本研究使用1︰50 000土壤矢量数据库实现栅格预测结果的输出。1︰50 000土壤矢量数据库是目前中国大区域尺度最为详尽的土壤数据库，由研究区第二次土壤普查各个县（市）图件资料GIS数字化而成。该数据库旱地图斑数目为17 024个，基本制图单元为土种，图例单元命名采用中国土壤发生分类系统。

1.3 样点数量设计

不同样点子集的创建采用ArcGIS 10.2地统计分析Create Subsets模块。首先从2008年可利用样点数据集中随机剔除197个，剩余884个，保证1980年和2008年样点总数一致。然后分别从两期的884个样点中随机抽取85%的样点（751个）作为训练数据集参与实际的空间插值过程，剩余的15%（133个）作为验证数据集用于插值预测精度的独立验证和不确定性评价。对751个训练数据集按5%的递减方式进行随机抽取，形成其余19组训练样点子集通过SPSS软件对各样点数据集的SOC进行最大值、最小值、标准差（SD）、均值、变异系数（CV）、偏度、峰度等描述性统计和正态分布检验。

1.4 空间预测方法及其精度评价

本研究采用当前土壤学使用最多和效果最稳健的普通克里金（Ordinary Kriging，OK）方法实现苏北旱地SOC含量的空间预测。空间预测结果输出为30 m×30 m的栅格，然后将每个土壤图斑所包含栅格像元（即SOC含量预测值）平均值提取至矢量土壤图图斑。

本研究选择常用的均方根误差（Root Mean Square Error，RMSE）和相关系数（r）评价不同样点数量下SOC空间预测的准确性。由于1980年和2008年SOC数据方差存在一定的波动性，RMSE和r不能直接用于不同时期预测精度比较。为了消除该影响，本研究也采用标准均方根误差（Normalized Root Mean Square Error，NRMSE）比较不同时期插值预测精度，其中：

式中，n为验证样点数量（本研究n＝133），Xmi和Xpi分别表示第i个验证样点中SOC含量的实测值和预测值。RMSE或NRMSE越小，r越大，预测精度越高。

1.5 经典统计学确定合理采样数目的计算公式

本研究所获取的最佳采样数量与经典统计学确定的合理采样数量进行对比，以分析其合理性。在目前大多数研究中，Cochran[14]针对与区域范围内纯随机采样而构造的计算公式通常被经典统计学用于最佳采样数量的确定：

式中，n为所需要的最佳采样数量，t为与显著水平对应的t式分布值，Std为样本标准差，d为样本平均值与允许误差的乘积。

2 结果与讨论

2.1 不同样点数量下1980年和2008年SOC含量统计学特征及空间变异结构

基于训练样点总集（751个样点）的SOC含量统计结果表明（表1），1980年研究区SOC含量均值为5.70 g·kg-1，而2008年含量均值明显提高，达到10.41 g·kg-1。两个时期其余19组样点子集SOC含量均值、变异系数与训练样点总集统计结果相似，表明随机抽取的不同样点数量均具有较好代表性，降低采样数量对研究区SOC基本统计信息影响较小，这与程道全等[15]和李凯等[16]的研究结果一致。进一步的Kolmogorov-Smirnov检验表明，1980年所有样点子集SOC含量全部符合正态分布，2008年样点子集SOC含量经对数转换后，也均符合正态分布，满足地统计插值的要求。

表1 不同样点数量下研究区SOC含量（g·kg-1）描述性统计结果及半方差模型拟合参数Table 1 Descriptive statistical and semivariance models fitting parameters of SOC contents (g·kg-1) in the study area relative to number of sampling sites

表征SOC空间变异的指示参数（块金值、基台值、结构性比例和变程）在不同时期或样点数量下具有较大差异，表明SOC空间变异结构随着样点数量以及时间的推移发生明显变化。1980年研究区SOC的结构性比例较小（3.43%～14.19%），其变异性受结构性因素影响在空间上强烈自相关；而2008年SOC结构性比例在14.29%～50.00%之间，绝大部分样点子集都超过了25%，土壤变异的空间相关性减弱为中等水平，这是因为1980—2008年研究区的耕作、施肥等土壤管理措施经历了较大改变，对SOC空间变异结构造成了一定影响。

2.2 不同样点数量下1980年和2008年SOC预测结果的空间分布特征

随着样点数量的减少，苏北旱地SOC的空间分布特征发生明显改变（图3和图4）。训练样点总集下，SOC空间分布既能反映整体格局又能体现变异细节，但随着土壤样点数量的降低，表征SOC局部变异的能力减弱，并逐渐改变整体空间分布格局，其中SOC含量高值和低值区域面积变化最为剧烈。造成不同数据集极值区域估算结果明显差异的原因可能是因为两端的极值具有较高空间异质性，样点数量减少导致表征土壤空间异质性的能力减弱。一方面，随着样点数量的降低，样点数据集的SOC含量最大值呈现减小趋势，1980年和2008年的变幅分别由751个样点的13.17和20.65 g·kg-1减小至38个样点的11.37和15.78 g·kg-1；而最小值呈增大趋势，1980年和2008年的变幅分别由751个样点的0.70和5.86 g·kg-1增大至38个样点的1.68和6.09 g·kg-1，表明样点数量减少会使得所探测到的SOC含量范围降低，进而影响其表征研究区的空间完整性。另一方面，巫振富等[17]研究指出，当最能凸显土壤空间分布格局的高、低值“关键区域”样点数量足够多时，才可保证预测结果最大程度地反映研究区土壤属性的实际空间分布特征，但样点数量的减少会导致“关键区域”能否满足样点需求成为一种不确定性事件。

续表

以训练样点总集空间预测结果的高、低值区为基准，划定研究区SOC空间分布的“关键区域”。经统计，751个训练样点总集下1980年和2008年SOC高值区的样点数量分别为147和40个，低值区分别为130和372个。当样点数量依次减少为638（85%）、563（75%）、451（65%）个时，位于“关键区域”的土壤样点数量仍然较多，其中1980年SOC高值区样点数量分别为121、116和97个，低值区样点数量分别为111、100和81个；2008年SOC高值区样点数量分别为33、28和29个，低值区样点数量分别为320、287和242个，尽管局部的一些细节信息逐渐被隐藏，但此时的SOC的整体分布格局和趋势仍然比较清晰，样点数量降低并未导致预测制图结果差异明显。然而，当样点数量进一步减少至113（15%）、38（5%）个时，位于SOC高值区的数量1980年分别仅有28和3个，2008年分别仅有8和2个；位于低值区的数量1980年分别仅有12和7个，2008年分别仅有48和17个，样点数量无法满足“关键区域”的制图需求，导致SOC在空间分布上的细节信息基本被平滑，呈现出单一的分布格局，不能很好地反映苏北旱地SOC的真实分布状况，特别是5%样点数量下空间结构出现失真畸变，表现的空间分布格局与训练样点总集相异。

从图2也可以看出，虽然2008年的采样点空间分布较不均匀，但该时期各个样点集的空间分布差异整体较小，这是由于其分布于“关键区域”的样点数较多，能更稳定地反映“关键区域”空间分布的真实情况，随着样点数量的降低，样点集统计特征的变异系数变化程度也明显小于1980年。这表明空间插值精度受采样数量和空间分布的共同影响，位于“关键区域”的样点数量越多，其估算结果越稳定。杨琳等[3]也认为，与单纯地增加样点数量相比，基于样点代表性的土壤采样在土壤空间制图成本和准确性上具有更为明显的优势。因此，在SOC高、低值等“关键区域”设置足够的样点数量，从而提高采样点所反映的空间异质性，是提高土壤空间预测效果的重要手段之一。

2.3 不同样点数量下1980年和2008年SOC空间预测精度

1980年和2008年研究区SOC含量空间预测精度随样点数量变化如图5所示。在1980年预测制图结果中，r和RMSE随样点数量变化的特征可划分为：100%～75%、70%～20%和15%～5%三个阶段。在第一阶段，r和RMSE变化幅度较小，预测精度并未随着采样点数量的减少而产生明显差异，说明当样点数量≥75%时，通过OK插值手段获得的SOC含量预测精度较高且已经达到相对稳定水平，研究区的最佳采样数目在75%（563个）左右，这可能是因为样点数量增加至一定程度以后，空间上的聚集度较高，导致信息达到一定的饱和水平。在这种情况下，即使采样点数量增加，也无法提供进一步的信息。在第二阶段，随着样点数量的减少，SOC含量的空间预测精度逐渐降低，r和RMSE变化幅度分别介于0.39～0.56和2.09～2.35 g·kg-1之间，原因在于预测精度受到样点代表性和密度显著影响。根据极值区域的样点统计，第二阶段的数量由216个减少至63个，SOC空间数据的代表性显著降低。此外，ArcGIS 10.2的空间统计分析表明，第二阶段样点之间的平均最近邻距离为4～9 km，明显大于前一阶段。在第三阶段，r和RMSE随着样点数量的减少急剧变化，预测精度较第二阶段降低得更快。半方差模型拟合结果显示，15%～5%样点数量下的拟合优度仅在0.08左右，SOC空间变异结构的推断存在较大误差，这可能是预测精度急剧降低的主要原因。然而，与1980年相比，使用20组相同样点梯度的训练子集数据对2008年研究区SOC含量预测时，预测精度随样点数量变化特征存在一定差异。2008年预测精度的变化趋势同样表现为先缓慢再剧烈下降，但阶段划分有所不同，主要分为：100%～70%、65%～20% 和15%～5%三个。在每个阶段，预测精度特征与1980年也不尽相同。第一阶段，预测精度并未随着样点数量的降低表现出显著变化，即70%的样点数时预测效果能达到相对稳定水平。可见，2008年研究区的最佳采样数目在526个左右，较1980年少5%，说明随着样点数量的增加，其预测精度更快趋于稳定。第二阶段，预测精度随样点数量降低呈现出明显的下降趋势，r和RMSE变化幅度分别处于0.42～0.63和2.11～2.56 g·kg-1之间。在该阶段内，2008年各个样点梯度间的预测精度差距明显大于1980年，表明样点数量的增加对其预测效果的提高更为显著。第三阶段，随着样点数量减少，样点梯度间预测精度差距进一步增大，即对样点数量的敏感性增强。但是，与1980年不同的是，2008年预测精度在15%样点数量下出现了一定回升，原因可能与样点数量稀少而引起的预测精度不确定性有关。研究表明，预测精度受到样点数量和空间分布的共同影响，随着样点数量的减少，空间分布位置对于预测效果的影响逐渐增大，甚至可以取代样点数量成为预测效果的主导性因素，但随机抽样忽略了样点的空间布局，使得样点是否对研究区SOC空间分布具有较好的全局代表性成为一种概率事件[17]。

从图5也可以看出，1980年和2008年苏北旱地SOC含量空间预测精度随样点数量变化表现出不同的响应特征，每个时期所需的样点数量有所差异。这主要是由于不同时期的环境、气候以及农田管理措施（施肥、灌溉和耕作等）不同，导致有机碳在土壤中的含量及其空间结构不同而产生响应特征差异。由于1980年研究区SOC空间自相关变程较小、局部变异较多，不仅需要更多的样点来推测空间分布，并且其预测精度的提高速率也较为缓慢和不稳定；相比而言，2008年研究区SOC空间自相关变程较大、局部变异较少，所以预测精度对于采样数量的敏感性较强，空间信息达到饱和状态所需样点数量也相对较少。从NRMSE来看，不同时期SOC含量的预测效果差距也十分明显，20组样点数量下1980年和2008年的NRMSE分别为0.34～0.43和0.20～0.25 g·kg-1，在样点数量相同的情况下2008年的空间预测精度优于1980年，主要的原因有以下3点：一是空间变异性对OK插值精度具有负面影响，该方法多适用于土壤属性变异较弱的环境区域[18]；二是由于不同时期样点的空间布局差异所致；三是2008年的土壤空间样点尽管分布不均，但位于关键区域的样点数量较多，其半方差模型拟合参数明显优于1980年，采样点具有更好的代表性。由此可见，研究区最优采样点数量及其预测精度不是固定的，而是由每个阶段土壤属性的空间变异特征以及样点的分布位置、空间布局共同决定。

根据经典统计学确定合理采样数目的计算公式（表2），在95%的置信区间和5%的误差条件下，1980年和2008年分别只需要261和100个采样点即可满足空间预测需求，而在90%的置信区间下，分别仅需20和8个采样点即可满足15%误差条件下的空间预测需求。本研究显示，1980年和2008年的最优采样数量应分别为563和526个，明显多于利用Cochran计算公式得到的最佳采样数量。Cochran在描述土壤元素空间变异时假定研究区土壤元素的含量变化是随机的，样本之间服从同一已知的概率分布并且相互独立，认为样品均值是反映研究区土壤特性的最优方式，决定土壤最佳采样数量的影响因素主要包括土壤样本的变异系数、允许误差和可接受的置信区间。然而，地统计研究表明土壤特性的变化并非完全随机，其在空间上具有一定的结构性和关联性，Cochran计算公式虽然能够显著减少采样数量，但只能概括SOC含量变化的宏观趋势而不能反映其局部变异特征，如果盲目采用Cochran公式确定采样数量，可能会造成超过预期的预测误差和不确定性。本研究采用随机抽样与地统计学的克里金插值相结合的方法确定合理采样数目，充分考虑了SOC的空间变异和结构特征，较Cochran将研究区作为一个均质的整体来确定最优采样数目的方式更具科学性。

表2 基于Cochran采样数量公式获取的苏北旱地SOC最佳采样数Table 2 Optimal number of sampling sites for prediction of SOC in the upland area in North Jiangsu based on the Cochran sampling quantity formula

3 结 论

本研究通过随机抽样方式设置20组不同数量的样点数据集，结合1︰50 000高精度土壤矢量数据库明确样点数量对1980年和2008年SOC普通克里金插值预测精度的影响，结果表明，1980年苏北旱地SOC预测值与实测值的相关系数和均方根误差分别在0.15～0.56和2.09～2.63 g·kg-1之间，样点数量约563个时可达到相对可靠的预测水平；2008年SOC预测值与实测值的相关系数和均方根误差分别在0.24～0.63和2.11～2.62 g·kg-1之间，空间信息达到稳定预测水平所需样点数量约526个。这表明相同区域范围内的最优采样数量及其预测精度不是固定的，是由每个阶段内土壤属性的空间变异特征以及样点的分布位置、空间布局共同决定。因此，未来的土壤采样设计需要根据研究区域土壤属性的时空变异特点制定合理的样点数量以及布局方式，以提高土壤采样的科学性和高效性。此外，有研究表明系统抽样能够较空间随机抽样更好地反映土壤属性的空间变异，未来的研究还应进一步考虑系统抽样方法对土壤属性变量空间预测结果的影响[19-20]。