黄淮海平原不同土地利用方式下土壤呼吸空间变异预测

李键,薛澄,2,杨扬,谢梦姣,彭正萍,2,王洋,2

(1.河北农业大学资源与环境科学学院/河北省农田生态环境重点实验室,河北保定 071000;2.河北农业大学国土资源学院,河北保定 071000)

土壤呼吸是全球土壤碳循环的重要组成部分,其空间异质性不仅能指示土壤微生物和理化性质的变化特征,还对估算区域碳收支起到重要作用。地统计学通过空间自相关作用揭示土壤属性变量的空间异质性,解释空间格局对生态过程与功能的影响[1]。此外,还能减少取样的不确定性和降低调查成本,其中普通克里金法已广泛用于土壤属性空间异质性研究[2],但其存在插值因子单一、不够精确的缺点[3]。协同克里金法则根据与主变量之间的逻辑关系和相关性引入多辅助变量,尤其当辅助变量与主变量高度相关时,协同克里金法在减少误差方面明显优于普通克里金法[4]。Golden等[5]对铅污染空间变异研究表明协同克里金法较普通克里金法均方根误差提高28.97%,王汉东等[6]在三峡区间面雨量插值研究表明协同克里金法较普通克里金法平均绝对误差提高19.28%,均方根误差提高10.25%。土壤呼吸具有较强的时空变异性[7],但应用地统计学方法预测土壤呼吸空间变异的研究较少,利用辅助变量提高土壤呼吸空间变异预测精度和有效性的研究仅在耕地进行了有限的探索[8],草地、林地等具有强大碳捕获潜力的土地利用类型中也亟待进一步开展更多实践。

基于辅助变量的协同克里金法可提高土壤呼吸预测精度。在不同土地利用方式下,由于植被生长及其主导的凋落物输入、土壤微生物种群差异造成土壤呼吸的不同,驱动耕地、草地、林地土壤呼吸空间变异的影响因子也存在差异。通常认为,土壤温度和水分是影响各土地利用方式下土壤呼吸季节变化的主导因子[9],但在空间异质性研究中,它们并不能完全解释土壤呼吸的空间变化。土壤有机碳可激发微生物的呼吸作用,其空间格局显著影响土壤呼吸的空间变异特征[10]。另有研究表明,土壤p H 值增加会引起土壤微生物类群改变,影响酶的组成和来源,对土壤呼吸产生显著改变[11]。此外,土壤N/P 变化,如磷添加措施会导致植被更易获取土壤磷素而减少植物根系向地下扩展,引起植物根系呼吸的降低[12-13]。总体而言,土壤呼吸空间异质性影响因子复杂多样,为充分解释土壤呼吸空间变异特征,辅助变量的筛选宜包含土壤温度、土壤水分、土壤p H 值及指示土壤养分特征的更多相关理化指标。

在本研究中,为探究不同土地利用方式下,土壤呼吸的空间格局以及确定哪些土壤环境指标能够在一定程度上辅助于土壤呼吸空间异质性的精确预测,选取黄淮海平原北部典型耕地、草地和林地,进行优化布点,采用普通克里金法和基于辅助变量的协同克里金法做插值分析,以探究:(1)不同土地利用方式下土壤呼吸的空间变异特征;(2)与普通克里金法相比,协同克里金法是否可提高土壤呼吸空间变异的预测精度。拟为今后研究不同土地利用方式下的土壤呼吸空间变异特征及精准预测区域碳收支提供数据支撑和理论依据。

1 研究区概况

耕地、草地、林地3种土地利用类型的研究样区均在河北省邢台市宁晋县贾家口镇域(114°46′—115°15′E,37°24°—37°48′N)范围内,位于黄淮海平原北部,太行山东麓的冲积平原上(图1)。研究区地势低平,属于暖温带大陆性气候区,年均气温约为13℃,年均降水量为450 mm。每年4 月中下旬到10 月上旬为无霜期,霜期160 d左右,研究样区土壤类型主要为潮土。

图1 研究区采样点分布示意图Fig.1 Distribution diagram of sampling points in the study area

2 材料与方法

2.1 采样点设计

在耕地、草地和林地分别设置50 m×50 m 的样方,基于6×6单元规则网格顶点设置规则网格样点36个(用×表示样点),辅以完全随机样点(用▲表示样点)和短距离样点(用+表示样点)的优化布点方案进行采样点设计(图1)。

耕地、草地和林地采样点总数分别为82 个、83个和79个,其中短距离样点均为12个。短距离样点的设置方法为随机选取1个已布设样点,设置在随机方向上与该点距离为δ的点作为短距离样点[14-15]。公式为:

式中:n为研究区采样点个数;L为每一个采样点与其最近采样点的距离,本研究中δ=4 m。

采样点密度越大、间距越小,短距离样点增加到10%时,可提升空间预测精度,有利于预测小尺度范围内空间异质性和可靠性[14-15]。采用S型多点混合采样法,取土壤样品深度为0—20 cm,将土样充分混匀,然后用四分法从混合后的土壤中留取500 g样品,除去砾石、根系、生物尸体等影响因素,经风干、研磨之后过筛装袋保存以待测定。

2.2 测定方法

2.2.1 土壤呼吸 于2018年9月夏玉米成熟期,放置LI8100测量室的PVC环(长10 cm,内径0.194 m),埋设深度为5 cm,以减少对表层土壤的扰动。测定前24 h剪除环内杂草,选取天气晴朗的上午,于9:00—11:00采用LI8100进行土壤呼吸速率测定,同时测定5 cm,10 cm土壤温度和5 cm土壤水分含量。

2.2.2 其他指标 土壤有效磷含量用钼锑抗比色法测定,土壤速效钾含量用乙酸铵浸提-火焰光度法测定,土壤全氮含量用凯氏定氮法测定,土壤p H 值用电极电位法测定,土壤有机质含量用重铬酸钾-外加热法测定,土壤有机碳含量根据有土壤有机质含量计算得来。转化公式为:

式中:SOC为土壤有机碳含量(g/kg);SOM 为土壤有机质含量(g/kg);1.724为van bemmelen因数。

2.3 半方差函数

2.3.1 普通克里金 半方差函数用于研究对象空间相关性,主要是在一维条件下,当采样点x在一维轴x轴上变化时,区域变量y(h)在点xi和(xi+h)处的值为y(xi)与y(xi+h)差的1/2,定义为区域变量y(x)在x轴上的变异函数,记为y(x,h)。公式为:

式中:h为变程;N(h)为以距离h为间距的实测对数;Z(xi)为样点xi处实测值。

2.3.2 协同克里金 协克里金法中使用的交叉半变异函数描述属性的空间变异性[3],传统的交叉半变异函数只有一个次级变量。公式为:

式中:y1为交叉半变异函数;h为变程;N为给定滞后距离区间(Xi+h)中Z1(Xi)和Z2(Xi)的个数。

基台值(C0+C)、块金值(C0)和变程(h)是半方差函数的参数,其中基台值表示变量在空间上的总变异大小,块金值表示随机因素引起的空间变异,变程反映空间变异的范围[16]。块金值与基台值的比值,即C0/(C0+C),为块金系数。当块金系数小于25%时,空间相关性最强;其值处于25%～75%,即为空间相关性中等;其值大于75%,则空间相关性较弱。

变异系数(CV)可反映土壤特性的空间变异程度[17]。当CV≤0.1时为弱变异,0.1＜CV＜1时为中等变异,CV≥1时为强变异。

2.4 空间插值方法

2.4.1 普通克里金法 克里金插值是一种基于变差函数理论和结构分析的有限区域内区域变量无偏最优估计方法。公式为:

式中:Z(x0)为位置x0的估计值;z(xi)为位置xi的已知值;λi为与数据相关的权重值。

2.4.2 协同克里金法 协同克里金法是普通克里金法的扩展应用,它利用了多种变量类型,将主变量的自相关和协同变量的互相关联结合在一起进行无偏最优估计。公式为:

式中:Z1(x0)为x0点的估计值;z1(xi)为初始变量xi的真实值;z2(xi)为同位二级变量xi的真实值;λi为与z1(xi)相关的权重;μi为与z2(xi)相关的权重;n1和n2分别为z1和z2的邻域个数。

2.5 精度评价

评价指标:平均绝对误差(MAE)、平均相对误差(MRE)、均方根误差(RMSE)和平均标准误差(MSE)。平均绝对误差确切地表示偏离真值的实际大小,平均相对误差反映模拟的可信程度,均方根误差反映模拟值的灵敏度和极值效应,平均标准误差表示样本均值与总体均值的离散程度。误差越小,效果越好。计算公式如下:

式中:zoi与zpi分别为土壤呼吸的实测值与预测值;n为样本数量。

2.6 数据处理

应用Microsoft Excel 2010计算均值、标准差,最大(小)值。SPSS 24.0做描述统计,K-S正态分布检验。Origin 2021制作相关性图。GS+(9.0)软件完成半方差函数模型拟合。ArcGIS 10.7 软件做空间插值、误差验证。

3 结果与分析

3.1 土壤指标统计特征

统计分析表明,耕地、草地和林地土壤呼吸速率的变化范围依次是1.08～3.54,1.75～6.35,0.50～3.77 μmol/(m3·s)(表1)。耕地、草地和林地土壤呼吸速率均呈现中等程度变异,变异系数依次为33%,24%,31%。

表1 土壤基本指标统计Table 1 Soil basic index statistics

土壤全氮、有效磷和速效钾含量在耕地,平均值分别为2.18 g/kg,39.59 mg/kg,319.94 mg/kg(表1)。大气温度、5 cm 和10 cm 土壤温度及土壤水分也表现为耕地最高,平均值依次为16.71℃,19.33℃,16.92℃,14%。土壤有机碳含量、p H 值、C/N 值在草地最高,平均值依次为10.99 g/kg,8.21,19.39。N/P和C/P在林地最高,平均值依次为0.09,1.55。3种土地利用方式下,大气温度、5 cm 和10 cm 土壤温度及土壤p H 值均表现为弱变异(仅耕地ST10和林地TEMP除外)。土壤水分、土壤有机碳、土壤养分及其相关化学计量学指标在耕地、草地和林地均呈现中等程度变异,其中变异最大的指标均为N/P,变异系数依次为50%,60%,67%。

3.2 土壤呼吸半方差和交叉变异函数

基于普通克里金法的半方差函数,在耕地和草地中为球面模型,林地中为指数模型(图2)。耕地和草地土壤呼吸均呈现强空间相关性,块金系数均小于3%,而林地土壤呼吸则呈现中等程度的空间相关性,块金系数为40%。土壤呼吸空间相关距离由远及近依次为林地(15.12 m)、耕地(3.83 m)和草地(2.47 m),这表明在区域范围内,林地土壤呼吸具有更远距离的空间相关性。

图2 土壤呼吸与协同变量的交叉变异函数Fig.2 Cross variogram of soil respiration and covariance variables

基于协同克里金法的研究表明耕地和草地土壤呼吸均为稳定模型,林地则是球状模型(表2)。土壤呼吸与辅助变量的交叉变异分析表示耕地、草地和林地的土壤呼吸与辅助变量之间都具有强空间相关性。仅耕地土壤中p H 值与土壤呼吸的空间相关性较弱(图3)。

表2 土壤呼吸插值模型误差Table 2 Error of soil respiration interpolation model

图3 土壤呼吸与因子相关性分析Fig.3 Correlation analysis of soil respiration and factors

3.3 土壤呼吸空间分布格局

基于普通克里金与协同克里金插值法预测的土壤呼吸空间分布格局大体保持一致,与普通克里金法相比,协同克里金法插值结果所表达的土壤呼吸空间格局在细节上更加准确、丰富,土壤呼吸高低峰值空间分布更加清晰(图4)。耕地土壤呼吸均呈现条带状,空间格局表现为四周高中部低。草地土壤呼吸空间分布均呈现不规则斑状,空间格局是以东北—西南为轴的对称分布。林地土壤呼吸空间分布均呈现明显的嵌套状,格局是自西向东逐渐增大。

图4 土壤呼吸速率空间插值Fig.4 Spatial interpolation of soil respiration rate

3.4 土壤呼吸空间预测精度评估

评估结果表明应用协同克里金法预测土壤呼吸空间变异的精度优于普通克里金法,其在耕地、草地和林地依次平均提升了3.27%,9.98%,7.86%(表2)。基于MAE,MRE,RMSE和MSE的预测精度评估表明,与普通克里金法相比,协同克里金法的土壤呼吸空间变异预测精度在耕地提高1.15%～6.34%,在草地提高9.10%～10.87%,在林地提高4.78%～12.09%(表2)。

4 讨论

4.1 不同土地利用方式下土壤呼吸异同

本研究发现,草地土壤呼吸速率较耕地和林地更高,这与前人的研究结果相一致[18-19]。Raich等[20]研究认为草地土壤呼吸速率大于耕地和林地,其指出土壤呼吸速率与不同土地利用方式下植被的差异密切相关,植被能通过影响土壤微生物、土壤有机质、土壤结构和根呼吸来影响土壤呼吸总量。有机物提供异养分解的基质[21],土壤有机碳含量作为环境背景值的重要参数,在一定程度上起到调节土壤呼吸作用。土壤有机碳储量越大的生态系统,土壤呼吸速率相对越高,因此研究样区中草地土壤呼吸速率比耕地和林地都更高。

4.2 辅助变量选取

选取辅助变量主要根据相关性和协同克里金模拟效果,辅助变量与土壤呼吸的相关性越大,预测效果越优[7]。耕地、草地和林地土壤呼吸速率均与温度(ST5,ST10,TEMP)和土壤水分呈现一定相关性(图3),即土壤呼吸速率随着温度或水分的升高而变大,这与前人研究结果相一致[17]。

草地和林地土壤呼吸的辅助变量除共有的土壤水分和5 cm 土壤温度外,林地的辅助变量还有气温,其与林地土壤呼吸呈显著正相关关系(p＜0.01)。耕地土壤呼吸的最优辅助变量中并未出现土壤温度和土壤水分指标,而是C/P,N/P和土壤p H 值,这可能与耕地长期施用有机肥、氮磷肥等措施有关,致使磷成为耕地土壤呼吸主导因子。孙宝玉等[13]研究表明施磷改变土壤N/P,土壤p H 值等理化指标,土壤中氮未受限制时,适量磷添加使植被更容易获取土壤磷素,相应地减少植物根系向地下的扩展,因此施磷引起的土壤N/P降低限制了植物根系呼吸。此外,施磷后凋落物碎屑、根系分泌物等有机碳输入数量和质量的变化,从而导致C/P 改变,进而影响土壤呼吸[22]。耕地中土壤p H 值提升,通过改变土壤酶空间构象来影响酶催化活性,显著影响土壤酶对碳底物的分解,来调控土壤呼吸[23]。因此,在耕地中,N/P,C/P和土壤p H 值为土壤呼吸空间变异的辅助变量。

4.3 预测精度评估

基于协同克里金法插值表达的土壤呼吸空间变异优于普通克里金法,可能归结于辅助变量对土壤呼吸空间变异的解释作用。与普通克里金法相比,协同克里金法的土壤呼吸预测在草地和林地中的MRE分别提升9.82%,12.09%,这与阳宽达等[24]研究结果相一致。胡丹桂等[25]采用RMSE 作为评价指标,在空气湿度空间插值研究表明协同克里金法预测精度提升7.66%～20.58%。本研究中,RMSE 提升4.19%～9.10%,与前人研究结果较为一致[25-26]。有研究表明,当主辅变量间的相关性超过0.45时,协同克里金法的预测精度明显优于普通克里金法[5],能够表达更为丰富的空间信息。而本研究中没有指标与土壤呼吸相关系数超过0.45,这可能是没有前人研究中精度提升那么高的重要原因,这值得进一步探究。

5 结论

(1)不同土地利用方式下土壤呼吸均呈现中等程度变异,平均速率由大到小依次为草地3.96μmol/(m3·s)、林地2.01μmol/(m3·s)和耕地1.97 μmol/(m3·s)。在空间自相关方面,空间相关距离呈现林地＞耕地＞草地,分别为15.12,3.83,2.47 m。

(2)本研究应用协同克里金法模拟土壤呼吸空间变异的最优辅助变量在耕地为C/P,N/P,p H 值,在草地为SM,p H 值,ST5,在林地为TEMP,SM,ST5。与普通克里金法预测相比,引入辅助变量的协同克里金法插值对耕地、草地、林地土壤呼吸空间变异预测精度的提升可达9.98%。因此,以土壤呼吸的影响因子作为辅助变量的协同克里金法可以更准确地预测区域土壤呼吸的空间变异特征。