土壤是地壳的风化层在气候、地形、时间、生物的共同作用下形成的非均质和变化的时空连续体^[1]，具有高度的空间变异性，即使在土壤类型和质地相同的区域内，同一时刻土壤特性(物理性质、养分含量等)在空间上也有明显差异，土壤的这种属性被称为土壤特性的空间变异性^{[2, 3, 4]}。土壤特性的空间变异是由土壤复杂的形成过程导致的，母质和气候在大尺度上决定了土壤性质的空间异质性，而在一个具体的生态系统内，土壤空间异质性除了受母质和气候的影响外，地形和生物的影响起主导作用，尤其是在小尺度上，地上地下的生物相互作用是导致土壤空间异质性的主要原因^[5]。不同尺度上的土壤空间异质性不但对土壤的结构与功能有重要的影响，而且对植被的空间格局也有重要的影响^[6]。在过去20年间，大量有关“异质性”、变异性”、“格局”和“尺度”方面的研究结果，迫使生态学家们放弃以往研究中所做的“理想化、均一、恒定的环境”假设^[7]。探讨特定尺度下的土壤空间异质性，对了解土壤性质的空间组成、空间构型和空间关联性，揭示空间格局对土壤生态过程的影响具有重要意义^[8]。然而，由于土壤物理性质本身的不均一性和空间变异性，以至环境和人为因素对土壤物理性质影响的规律或变化趋势，人们仍不是十分清楚^[9]。

蜀南竹海是我国唯一以大面积竹景为特色的国家级风景名胜区，也是我国毛竹富集区之一。20 世纪80 年代末，蜀南竹海景区毛竹林成片死亡，是该区林地土壤退化的明显标志，是影响该区毛竹产业发展的瓶颈^{[10, 11]}，恢复与提高毛竹林生产力已成为当务之急。现有毛竹林土壤的研究，大多集中在土壤养分动态、渗透性能以及不同经营类型毛竹林土壤物理、化学和生物性质(含土壤呼吸)等方面^{[12, 13, 14, 15, 16]}，至于毛竹林土壤物理性质的研究，侧重于剖面的垂直分异或者不同林型条件下的变化，缺乏土壤物理性质空间异质性的深入研究。本文以蜀南竹海毛竹林土壤为研究对象，应用地统计学方法，研究其物理性质空间异质性，以期为维护与提高毛竹林生产力，增强生态服务功能提供理论和试验依据。

蜀南竹海风景名胜区(26°33′17″—28°26′46″N，104°5′11″—105°4′54″E)位于四川省宜宾市长宁、江安两县连接地带，地处四川盆地南缘向云贵高原过渡的连天山余脉。景区地势中部最高，东部次之，西部最低，海拔高程在260—1000 m之间。总面积达120 km²，核心景区42.228 km²。研究区属亚热带湿润季风气候，年降水量1200—2000 mm，年平均气温14.5—18°C。景区内主要的植物为竹类，覆盖率大于87%，其中毛竹(Phyllostachys heterocycla cv pubescens)占91%，其他还有凤尾竹、琴丝竹、方竹、慈竹、黄竹、水竹、黑竹、斑竹等，亩均立竹170根以上，该区毛竹林平均年龄为6a左右。出露岩层属白垩系夹关组紫红色厚层长石砂岩、粉砂岩、页岩。土壤类型为山地黄壤和紫色土，质地大多为砂壤，透水性能好，土壤pH约5。

以核心景区1∶10000地形图为工作底图，按照公里网格(1000m×1000m)设置34个样点(图 1)。用GPS导航，到实地确定样点位置。若图上标明的位置在野外无法到达，或为水域、农田等时，就近另选样点，并在图上标明准确位置，重新用GPS确定其经纬度。在每一样点附近挖3个临时剖面，作为采样点，记录各采样点坡向、海拔、土层厚度等(表 1)。由于土壤的生物活性主要集中在表层几cm—30 cm^[17]，因此，每个采样点按照0—20 cm、20—40 cm土层分别用环刀取样，供分析土壤的水分物理性质。同时另取50 g土样放入铝盒中，分析土壤水分含量，再取200 g土样装入塑料袋中，供土壤粘粒含量分析。

图 1 样点分布图及研究区边界 Fig.1 Scattergram of border of study area and sampling sites

图选项

土壤物理性质是衡量土壤肥力和林地环境质量高低的重要参数。本文选取土壤粘粒、容重、毛管孔隙度、非毛管孔隙度、毛管持水量、最大持水量和最小持水量等7个因子，作为度量土壤物理性质的参数。土壤水分含量采用“烘干法”测定，容重、毛管持水量和孔隙度采用“环刀法”测定，土壤粘粒含量采用“吸管法”测定^[18]。每个样点附近3个临时剖面土壤物理性质测定值取其均值，为该样点土壤的物理性质。

(1)半变异函数

半变异函数(semivariograms)是地统计学中用以探索土壤空间变异性的最关键的函数，也称为半方差函数，反映了土壤性质在不同距离上观察值之间的变化。半方差函数是一个连续的函数，描述了土壤性质空间变异的连续性^[19]。在一维条件下，当空间点x在一维x轴上变化时，区域变量Z(x)在点x和(x+h)处的值Z(x)与Z(x+h)差的方差的一半，定义为区域化变量Z(x)在x轴方向上的变异函数，记为γ(h)，即

点击浏览公式

式中，γ(h)是半方差函数，h为样本间距又称为步长或位差(Lag)，N(h)表示的是间距为矢量h的所有观测点对的数目。Z(x_i)、Z(x_i+h)表示Z(x)在x_i和x_i+h位置上的实测值。图 2为变异函数曲线示意图。

图 2 球面模型半方差函数曲线示意图 Fig.2 Semivariogram of spherical model

图选项

图中，C₀称为块金值，它表示h很小时，两点间取值的变化，反映了区域化变量Z(x)内部随机性的可能程度；a为变程，它的大小反映了研究对象中某一区域化变量的变化程度，也反映了影响范围。基台值由拱高C和块金值C₀之和组成，反映某区域化变量在研究范围内变异的强度。通常使用块金值与基台值之比C₀/(C₀+C)或拱高与基台值之比C/(C₀+C)，反映随机部分或自相关部分占总空间异质性的大小。根据Cambardella等^[20]研究结果，当C₀/(C₀+C) < 25%时，表示土壤物理性质有很强的空间依赖程度，25%≤C₀/(C₀+C)≤75%时，表示土壤物理性质有中等程度的空间依赖程度，C₀/(C₀+C)>75%时，相应的土壤物理性质被认为有很弱的空间依赖。决定系数表示变异函数与模型拟合的程度，决定系数越大，模型拟合度越高^[21]。

(2)空间插值方法(克立格，Kriging法)

空间插值是使用地统计学方法，以采样间隔比较大的土壤质量指标数据，对未采样点进行插值，从而得到面状的数据，以便生成土壤质量等级图。由此可知，土壤质量精确管理的一个重要基础，是较理想的空间插值结果。用于插值分析的方法很多，克立格(Kriging)插值法是目前应用比较广泛的一种^{[22, 23]}。克立格插值实际上是一个进行局部估计的加权平均值，其表达式为

点击浏览公式

式中，Z(x_o)是未采样点x_o处的内插估测值，Z(X_i)表示的是x_o附近的若干采样点的实测值。

采用经典统计方法，计算土壤物理性质的平均值、最大值、最小值、标准差、变异系数、偏度、峰度、W检验等。应用地统计方法对土壤物理性质进行分析，建立变异函数理论模型，并通过块金值、基台值、变程、结构比等参数，分析土壤物理性质的空间变异规律，再运用克立格空间局部插值估计法(Kriging)，对研究区内不同位置的观测指标进行插值估计，揭示空间分布格局。

基本统计分析利用SPSS 17.0 完成，地统计学分析用GS+ Version9完成，克立格空间局部插值估计用arcGIS9.3空间分析模块中的普通克立格法(Ordinary Kriging)完成。

经典统计分析表明(表 2)，研究区0—20 cm土壤粘粒、非毛管孔隙度、毛管孔隙度、毛管持水量、最大持水量和最小持水量均大于20—40 cm土层，惟土壤容重正好相反。说明研究区土壤随土层厚度增加，土壤粘粒、非毛管孔隙度、毛管孔隙度、毛管持水量、最大持水量和最小持水量减小，而土壤容重则随土层加深而增大。按照反映离散程度变异系数(CV)的大小，变异系数可以粗略分级为弱变异CV < 10%，中等变异10%≤CV≤100%，强变异CV>100%^[24]。因此，从表 1可知，0—20 cm和20—40 cm土壤变异系数变化范围分别为15.17—40.87%和14.78—57.09%，均为中等程度的变异；并且0—20 cm土壤毛管持水量、最大持水量、最小持水量、毛管孔隙度、容重的变异均大于20—40 cm的变异，而土壤粘粒、非毛管孔隙度变异均小于20—40 cm的变异。这与土层越深，土壤石砾含量增多以及毛竹根系分布状况等有很大关系。

经典统计学中的变异系数仅仅从统计学的角度描述了这些土壤物理性质的变化程度，但不能反映这些物理性质的空间分布格局与空间异质性尺度。表 3 是土壤物理性质变异函数的理论模型及有关参数。从表 3可以看出，物理性质空间变异尺度即变程为404.9—11887.1 m。其中，0—20 cm土壤粘粒的变程最大(11887.1 m)，0—20 cm土壤最小持水量的变程最小(404.9 m)；0—20 cm土壤粘粒、非毛管孔隙的空间变异尺度大于20—40 cm土层，而毛管孔隙度、容重、毛管持水量、最大持水量和最小持水量的空间变异尺度小于20—40 cm土层。因为两个土层土壤最大持水量(636.0、710.7 m)和最小持水量(404.9、534.1 m)的变程小于采样间距1000 m，所以不符合内插条件，难于进行空间插值计算，而其他物性质变程均大于采样间距1000 m，符合内插条件要求，可以进行空间插值计算。

从表 3还可以看出，0—20 cm土壤物理性质的基台值(C+C₀)均大于20—40 cm土层，说明0—20 cm土壤物理性质总的空间变异程度大于20—40 cm土层。另一方面，0—20 cm非毛管孔隙度和20—40 cm毛管孔隙度的C₀/(C₀+C)较大，在25—75%之间，属于中等程度的空间依赖程度，而两个土层的其他物理性质均小于25%，属于很强的空间依赖程度，特别是两个土层的土壤容重和20—40 cm土壤粘粒空间自相关性较大(基底效应小于1%)。说明研究区虽然人为等随机干扰较严重，但大多数土壤物理性质稳定性较强，具有较好的空间结构性。

根据以上分析结果可知，研究区土壤粘粒、容重、毛管持水量、非毛管孔隙度和毛管孔隙度等5个物理性质符合空间内插条件，于是按照其变异函数模型进行克里格插值，分别绘制0—20cm和20—40 cm空间分布图(图 3,图 4)

图 3 毛竹林土壤0—20 cm土壤物理性质空间分布格局 Fig.3 Spatial distribution of soil physical factors in the 0—20 cm on bamboo forest

图选项

图 4 毛竹林土壤20—40 cm土壤物理因子空间分布格局 Fig.4 Spatial distribution of soil physical factors in the 20—40 cm on bamboo forest

图选项

从图 3、图 4看出，土壤粘粒含量在0—20 cm土层分布格局较20—40 cm明显，并且由东南角粘粒较少的斑块呈放射状向东北和西北递增，并在东北角和西北角达到最大值。图 3、图 4中，土壤容重较小的斑块都出现在中部地区，而在北角容重都偏大。其中，0—20 cm土壤容重空间分布呈相互交替的不规则点状，20—40 cm土壤容重空间分布格局较0—20 cm明显，斑块之间过渡明显，但两个土层土壤容重均在东北角和西北角较大。图 3-c、图 4-c表明，两个土层土壤毛管持水量空间分布格局较为相似，呈较大的长块状贯穿于南北之间并向东北角延伸，较小的分布在西北角，这与研究区地势中部高，b东西部较低的变化相似，并与赵超^[25]在不同海拔毛竹林土壤特征的研究结果相一致。由图 3、图 4可看出，两个土层土壤非毛管孔隙度空间分布格局均较明显，两个土层之间也较相似。非毛管孔隙度较小的两个斑块镶嵌在较大的基质上，斑块之间有明显的过渡现象。但不同土层之间也有一定的差异，0—20cm土层非毛管孔隙斑块组成更大，过渡现象更明显，非毛管孔隙度较少的分布在中南和中西部，这与0—20cm土壤容重分布格局正好相反。说明它们之间存在直接或间接制约作用。图 3、图 4所示，两个土层土壤毛管孔隙度空间分布不明显。 其中，0—20 cm土壤毛管孔隙度最大的斑块分布在北部和南部，最小的分布在西部；20—40 cm土壤毛管孔隙度较大的斑块分布在中部，整体上由中部向东西两个方向递减。

从表 4可看出，土壤毛管孔隙的空间变异与毛管持水量、非毛管孔隙度的空间变异有密切关系，并且达到极显著水平；土壤容重与粘粒的空间变异有密切关系，并且在20—40 cm土层达到极显著水平；20—40 cm土壤容重与非毛管孔隙的空间变异有密切关系。说明随着土层厚度增加，各物理性质之间的相互作用越突出。这主要是因为深层土壤受到外界干扰影响较少，土壤性质之间的相互作用占主导作用。

近年来，对土壤物理性质空间异质性的研究逐渐受到重视，但由于不同的生态系统和不同的采样尺度，所得土壤物理性质空间异质性的研究结果有很大的差异。李岩^[26]等人在阔叶红松林林隙土壤物理性质微环境(采样间距为2—22.6 m)异质性研究中发现，土壤毛管持水量的变程在0—20 cm为4.39 m，空间依赖程度为15.10%；20—40 cm变程为1.83 m，空间依赖程度为21.90%。Iqbal等^[27]研究了农田土壤物理性质的空间异质性，发现表层土壤容重的变程为106m，空间依赖程度为29%。与本研究相比，李岩和Iqbal等的采样尺度和最小间隔距离均小于本研究(1000m)，因此得出的变程也小于本研究，而孙玉冰等^[28]的研究尺度(研究样地为1141km²，最小采样间隔为8 km)远远大于本研究尺度，所得变程(容重变程11.90 km)也远大于本研究。这是因为在不同的尺度上，影响土壤物理性质空间分布格局的因素和生态过程不同，即微小尺度上土壤微生物活动为主要影响因素；在cm到m的尺度上，土壤动物的活动为主要影响因素；而在更大尺度上，地形、土壤母质、植被等对土壤物理性质的空间结构影响较大^[29]。因此，土壤物理性质的空间结构特征随采样间距大小的变化而变化，在特定采样尺度下便显示相应的变化规律。

在本研究中，土壤粘粒、容重、非毛管孔隙度、毛管孔隙度、毛管持水量等5个物理性质的变程为1252.3—11887.1m，均大于1000m，因此采用公里网格布点是合理的。土壤最大持水量和最小持水量的空间变异尺度较小(404.9—710.7 m)。这是因为土壤最大持水量受环境因子影响较大，特别是坡向、坡度、地貌、海拔和土层厚度等因子影响^[30]，土壤最小持水量与土壤结构、质地、有机质含量以及土地经营状况有关，不同土壤的最小持水量不同，即使是同一种土壤，由于经营方式不同，最小持水量亦不同^[31]。因此，在以后类似土壤水分物理性质研究中，采用公里网格布点不能较好反映他们的空间分布结构，应适当缩小采样间距，不能超过500 m，以300—400 m为宜。这是因为如果采样间距过大，难于精确反映土壤物理性质的空间结构性，而采样间距过小，则造成人力、物力的浪费。土壤调查与研究中，布点抽样是一项基础性的工作，应用地统计学方法研究土壤物理性质的空间异质性，如何在保证精度的前提下确定适宜的采样间距，需要进一步的研究。

土壤物理性质空间变异与许多因素有关，是母质、地形、气候等结构性因素和土壤生物、人为干扰等随机因素共同作用的结果。研究区土壤的物理性质，除土壤容重外，其数值随土层厚度增加而减小，土壤物理性质的变异系数，总体上也是0—20cm土层大于20—40cm土层，但不同土层的土壤物理性质都存在空间自相关性；另外，从土壤物理性质的空间分布格局可以看出，研究区中部土壤的物理性能优于东部和西部。这一方面与研究区地势中部较高，东、西部较低有关(与赵超等人的研究结果一致)，另一方面与风景区活动中心在东、西部，游客满载，加上近年来工程建设，频繁的采竹、挖笋等人为干扰有关。

森林对土壤一些性质的影响，首先表现在土壤物理性质的改变上。这是因为林木着生土壤时，首先在土壤表层发展强大的根系，给通气和透水创造良好条件，改变持水量，改善土壤的结构等土壤水分物理性质^[32]。因此，研究区一要科学管理，尽量减少人为干扰，二要优化林分结构，改善土壤物理性质，增强土壤的生物自肥作用，维护与提高毛竹林生产力。