邓佳楠，刘 帆，刘泽彬，郭建斌，王彦辉
（1.北京林业大学水土保持学院，北京 100083；
2.中国林业科学研究院森林生态环境与自然保护研究所/国家林业和草原局森林生态环境重点实验室，北京 100091）
土壤密度是土壤的基本物理性质之一，对坡面土壤水分运动有重要影响[1]。受气候（如降雨）、土壤质地、地形和植被类型等差异的影响，土壤密度在坡面上有很大的空间异质性。准确探明坡面土壤密度空间变异及其影响因素对于深入理解坡面产流机制和土壤水分的坡面分配有重要意义[2]。
土壤密度空间变异研究涉及到了区域、流域、坡面和样地等不同的空间尺度[3-7]以及草地、灌丛和森林等不同的植被类型[8-11]，但关于坡面土壤密度空间变异的研究，主要集中在草地[8]和灌丛[10]。现有的森林坡面研究多关注土壤密度的坡位差异[12]，但通过几个坡位的简单比较还不能完整揭示坡面微地形、立地环境等各差异下的土壤密度空间变异特征，这是因为即使在特定坡位的土壤密度也存有较大空间变异[11]。因此，森林坡面土壤密度空间变异研究还需细化。影响土壤密度空间变异的因素在不同空间尺度存在差异，在坡面上，地形（如坡度）、植被类型和土壤质地（结构）等是影响坡面土壤密度空间差异的重要因素[13]。但在不同地区、植被类型下，其主导因素也有差别，如影响土壤密度空间变异的主导因素在河北坝上地区草地坡面为地形[8]，在四川西缘山地坡面是植被[14]，在喀斯特坡面是地质[10]。作为黄土高原的产流核心区-土石山林区，土壤密度异质性是影响坡面产流空间变化的重要因素，但坡面土壤密度空间变异的影响因素未有探索。另外，多数研究关注表层或特定深度土壤密度的空间变异[10,15-16]，但受不同土层石砾含量、根系生物量等差异影响[17-18]，不同土层土壤密度空间变异可能存有差异，但目前关于土壤密度空间变异随土层深度变化还缺乏足够探索。
为克服坡面土壤密度空间变异的影响，往往需进行高密度的取样来确保坡面土壤密度的评估精度，但费时、费力[19]。若能量化土壤密度合理取样数及与空间变异的关系，将显著提高坡面土壤取样的效率和评估精度。目前有关土壤密度合理取样数的研究较少，刘帆等[11]虽量化了华北落叶松林表层土壤密度的合理取样数，但这仅是样地尺度的探索，且未探讨土壤密度空间变异与其合理取样数的关系。
华北落叶松（Larix gmelinii var.principisrupprechtii （Mayr） Pilg.）是黄土高原土石山区主要的造林树种，对当地的水源涵养、保持水土有重要作用。理清华北落叶松林坡面土壤密度空间变异及合理取样数对于准确理解华北落叶松林坡面产流空间差异的影响机制具有重要价值。为此，本研究以六盘山半湿润区的华北落叶松林坡面为研究对象，结合经典统计学和地统计学方法，明确华北落叶松林坡面0～100 cm土层土壤密度的空间异质性及影响因素，并采用Monte Carlo模拟重抽样方法，量化各层土壤密度的合理取样数。
1.1 研究区概况
研究区位于黄土高原土石山区的六盘山香水河小流域（106°12′10.6″～106°16′30.5″ E，35°27 ′22.5 ″～35°33 ′29.7 ″ N），流域面积为43.7 km2，海拔2 010～2 942 m，土壤类型以山地灰褐土为主，石砾含量较高；
气候为暖温带半湿润气候，年均气温6.0 ℃，年均降水量618 mm。小流域内主要以人工林和天然次生林为主，其中华北落叶松是主要的人工林树种，占人工林总面积的90%。
在小流域内选取了1个东南坡向的典型华北落叶松人工纯林坡面，坡面上华北落叶松均为同龄（39 a，2020）林。海拔范围2 259～2 478 m，坡长480.6 m，平均坡度为27.8°，土层深度为1.0～1.2 m左右，土壤总孔隙度55.2%，石砾含量9.5%。林分密度为822株·hm−2，郁闭度0.73，平均树高17.1 m，平均胸径19.8 cm，生长季平均叶面积指数2.90。林下灌木零星分布，以秦岭小檗（Berberis circumserrata （Schneid.） Schneid.）、蒙古荚蒾（Viburnum mongolicum （Pall.） Rehd.）和西北栒子（Cotoneaster zabelii Schneid.）为主。草本覆盖度在40%左右，以东方草莓（Fragaria orientalis Lozinsk.）和点叶薹草（Carex hancokiana Maxim.）为主。
1.2 土壤密度、孔隙度和石砾含量测定
2020年9月，参考刘春利等[20]关于坡面土壤水文物理性质空间变异的布点方法，沿坡面设立3条间距为15 m的样线，在每条样线自下而上每隔30 m设置1个取样点，共48个。在每个样点按 0～20、20～40、40～60、60～80和80～100 cm 5个土层，用体积为200 cm3的环刀在各层的中部竖直向下取原状土，并记录取样点坐标（以最左侧样线下方第一个取样点为原点）。将样品带回室内后，测定土壤密度、总孔隙度等，用排水法测定石砾含量。
1.3 空间变异分析
变异系数（CV）是经典统计学中分析空间变异程度的衡量参数。根据变异程度强弱，CV可划分为3个等级：弱变异（CV≤0.1）、中等变异（0.1 半方差函数是地统计学中分析空间结构的主要工具。半方差函数的主要参数为块金值（C0）、基台值（C0+C）和变程（a）。块金值C0表征随机变异程度； 利用Monte Carlo随机抽样方法，确定在一定误差下坡面土壤密度的最少取样数。首先，分别将 0～20、20～40、40～60、60～80和80～100 cm土层的土壤密度测值作为Monte Carlo随机抽样的初始值； 半方差函数分析在GS+9.0软件中进行，土壤密度Kriging插值采用Arc GIS 10.3软件，Monte Carlo模拟抽样在R 3.6.1软件中进行。 由表1可知，随土层加深，土壤密度逐渐增加，0～20、20～40、40～60、60～80和80～100 cm土层土壤密度平均值（变化范围）分别为0.88（0.73～1.08）、1.04（0.78～1.55）、1.23（0.91～1.66）、1.41（1.02～1.78）和1.49（1.09～1.73）g·cm−3，变异系数为0.09、0.16、0.17、0.15和0.11，土壤密度在0～20 cm土层为弱变异，在其它土层为中等变异。由图1可知，土壤密度空间变异系数与石砾含量的空间变异系数呈显著的（P＜0.05）线性相关，表明土壤密度空间变异随土层深度的变化主要取决于石砾含量的空间异质性。 图1 土壤密度变异系数与孔隙度和石砾含量变异系数的关系Fig.1 Relationships between the coefficient of variation （CV） of soil bulk density and the CVs of soil total porosity and rock fragment 表1 坡面土壤密度的统计特征Table 1 Descriptive statistics characteristics of soil bulk density on the slope 2.2.1 土壤密度的半方差函数分析 由表2可知，土壤密度的半方差函数在0～20和20～40 cm土层为高斯模型，在40～60和80～100 cm土层为球状模型，在60～80 cm土层为指数模型。各层土壤密度块金值均为正值； 表2 坡面不同土层土壤密度的半方差函数理论模型及其相关参数Table 2 Semi-variogram theoretical models of soil bulk density in different soil layers and their related parameters 2.2.2 坡面土壤密度的空间分布 由图2可知，坡面各层土壤密度有明显的空间异质性，自坡顶至坡底，土壤密度总体上呈“低～高～低”的变化。各层土壤密度的最大值基本在坡中，最小值在坡下（0～20 cm除外，在坡上）。0～20 cm土层土壤密度跨度范围较小，斑块较大，空间分布较均匀； 图2 坡面 0～20（a）、20～40（b）、40～60（c）、60～80（d）和80～100（e）cm 土层土壤密度的空间分布Fig.2 Spatial distribution of soil bulk density in 0～20（a）,20～40（b）,40～60（c）,60～80（d） and 80～100（e） cm soil layer on the slope 基于Monte Carlo模拟重抽样的方法确定了在95%的置信区间下，若控制坡面土壤密度在10%误差内的最少取样数。由图3可知，随取样数增加，土壤密度的平均值和置信区间逐渐向内收缩。在95%置信区间下，当取样数分别为4、10、10、8和5个时，得到的0～20、20～40、40～60、60～80和80～100 cm土层土壤密度可控制在10%误差内。由图4可知，土壤密度的最少取样数与其变异系数呈极显著的线性相关（P＜0.01），即随土壤密度空间变异系数的增加，最少取样数逐渐增加。 图3 Monte Carlo 模拟的 0～20（a）、20～40（b）、40～60（c）、60～80（d）和80～100（e）cm土层土壤密度平均值和置信区间随取样数的变化Fig.3 Variations of the average and confidence interval of soil bulk density in 0～20（a）,20～40（b）,40～60（c）,60～80（d） and 80～100（e） cm soil layer simulated by Monte Carlo as the sample size 图4 坡面土壤密度最小取样数随其变异系数的变化Fig.4 Variation in the minimum sample size of soil bulk density with its coefficient of variation 研究坡面100 cm土层土壤密度均值为1.21 g·cm−3，低于全国100 cm土层土壤密度均值（1.32±0.21 g·cm−3）[3]，这是因森林根系的作用增加了土壤的疏松性和通气性，使土壤密度降低[22]； 各层土壤密度的块金值均较小，表明随机因素对土壤密度空间变异影响较小。本研究坡面土壤密度的基台值与样地尺度[11]的基台值大致相当，表明在中小尺度上土壤密度总的空间变异相近。土壤密度的总空间异质性在0～20 cm土层最小（基台值最小），在40～60 cm土层最大（基台值最大），这与经典统计学的研究结果（变异系数）相吻合。各层土壤密度块金比在0.25～0.75间，表明各层土壤密度均为中等空间自相关，这与周卓丽等[8]在草地坡面研究得到的不同土层土壤密度的空间自相关性存在明显差异的结论有所不同，可能是与不同植被坡面的地形（如坡度）、立地（孔隙度、石砾含量、机械组成等）、生物因素（如根系生物量等）等空间变异不同有关[25-26]。随土层加深，土壤密度的变程先减后增，可能是因为土层中部受根系、石砾含量等因素的综合作用更强，导致土壤密度的空间自相关范围小。土壤密度沿坡空间分布在不同植被类型坡面存有较大差异，如杨宽等[27]得到福建高山毛竹林上坡位和下坡位土壤密度高于中坡位； 本研究发现，土壤密度最小取样数与其变异系数呈显著相关，表明土壤密度最小取样数主要取决于其空间变异程度。刘帆等[11]研究表明，当林地（样地）土壤密度空间变异系数为0.11时，最少取样数为4个； 在本研究的华北落叶松林坡面上，在0～20 cm土层土壤密度为弱变异，在其他土层均为中等变异，不同土层土壤密度空间变异的差异与石砾含量的空间变异密切相关。土壤密度的最少取样量取决于其空间变异程度，当坡面各层土壤密度空间变异系数在0.09～0.17时，若要控制100 cm土层土壤密度在10%误差内（95%的置信区间下），至少需要布设10个以上的取样点。
基台值（C0+C）为半变异函数的最大值，该值越大，表明总的空间异质性越高[10]；
块基比[C0/（C0+C）]用于判断变量空间自相关程度，块基比<0.25，表明变量呈强空间自相关性；
0.25～0.75，表明变量呈中等空间自相关性；
> 0.75，表明变量呈弱空间自相关性[21]。变程（a）是变量空间自相关范围大小的参数，在变程范围内，两点空间距离越小，空间自相关性越大；
当两点空间距离大于变程后，两点的空间自相关性不存在。1.4 合理取样数确定
然后对每个土层48个土壤密度测值进行不放回抽样（抽样数n=2～47），为了保证充足的随机样本，每个抽样数下（如n=2）重复1 000次，计算不同抽样数下土壤密度均值；
最后，计算并绘制不同抽样数土壤密度的平均值和置信区间（95%）。根据土壤密度均值及其置信区间随抽样数的变化，确定在10%误差内各土层土壤密度的最少取样数。1.5 数据分析
2.1 坡面土壤密度基本统计特征
2.2 坡面土壤密度的空间结构特征
土壤密度基台值沿土层深度的变化范围0.011～0.060，在0～20 cm土层最小，在40～60 cm土层最大。0～20、20～40、40～60、60～80和80～100 cm土层的块基比分别为0.295、0.348、0.300、0.495和0.498，表明各层土壤密度均属中等空间自相关性。各层土壤密度变程分别为257.70、126.20、179.30、215.42和294.50 m，表明随土层加深，空间自相关范围先减小后增大，在80～100 cm土层最大，在20～40 cm土层最小。
40～60 cm土层土壤密度跨度范围较大，条带密集，空间分布格局最复杂。随土层加深，土壤密度空间分布从斑块状逐渐变为条带状。
2.3 坡面土壤密度的最少取样数
另外，本研究区为土石山区，石砾较多，土壤孔隙度较高，导致土壤密度偏低。通常，土壤密度随土层深度增加而增加，本研究也得到同样的规律，其原因主要有两方面：一是土壤有机质的影响，研究表明[23]，土壤密度大小与土壤有机质含量显著负相关，土层越深，土壤有机质含量越低，土壤密度则会越大；
二是土体重量压缩的结果，土层越深，土体越紧实，土壤密度越大。这与周武等[24]得出的研究结果一致。土壤密度的变异系数在0～20 cm土层为弱变异，其他土层为中等变异，可能与不同土层石砾含量的空间异质性有关（图1），0～20 cm土层土壤主要来自于腐殖质和土壤有机质，石砾含量较少，土壤较为均匀，空间变异小；
其他土层石砾含量较多，且分布不均，空间变异大。
高雪松等[14]得到四川西缘山地坡面土壤密度表现为上坡位>中坡位>下坡位。本研究中，土壤密度从坡顶至坡脚总体呈“低～高～低”的变化，在坡中最大，这种差异可能与不同坡面土壤质地、地形、植被类型和生长特征、人为干扰等差异有关。
本研究中，当坡面80～100 cm土层的土壤密度空间变异系数为0.11时，最少取样数为5个。这表明在同等空间变异下，空间尺度越大，所需的最少取样数可能越多，但该结论未来还需更多数据进行验证。由于各层土壤密度空间变异不同，导致不同土层土壤密度最少取样数存在较大差异，在进行野外采样时，应基于取样深度及取样深度内土壤密度的空间变异程度制定合理的取样策略。本研究坡面，为保证100 cm土层土壤密度的评估精度，至少需布设10个以上的取样点。