青藏高原是我国主要牧区之一, 草地类型丰富, 其面积占我国草地面积的42%^[1]。近年来, 由于环境变化以及人类生产活动的影响, 青藏高原草地退化日趋严重, 生产力急剧下降^[2]。过度放牧是导致草地退化的主要因素之一^[3], 不仅抑制植被生长, 并且降低土壤肥力^[4-7]。因此, 必须采取有效的生态措施来恢复退化的草原生态系统, 如围封禁牧、禁牧+补植和合理放牧^[8]等措施, 目前, 关于不同利用方式对退化草地生态系统恢复的影响研究越来越多, 但并没有得出一致的结论。仲波^[9]对若尔盖草地土壤微生物的研究中发现, 适度放牧提高了酶活性以及土壤碳氮矿化速率, 促进退化草地的恢复。刘艳萍^[10]在研究不同恢复措施对退化草地土壤和植被的影响中发现, 禁牧显著改变群落结构和增加土壤养分。

土壤粒径分布属于土壤物理性质, 可以表征土壤结构, 反映土壤的发育状况, 对土壤水、肥等物理特性具有一定的影响。土壤是由形状与大小各异的土壤颗粒组成的多孔介质^[11], 具有一定的分形特征。分形理论作为一种揭示物质局部结构和形态的工具在研究土壤结构中得到了广泛运用。例如, 孙梅等^[12]研究长期不同施肥条件下红壤粒径分维, 表明施用有机肥可以促进土壤细化改善土壤粒径分布格局。王德等^[13]对黄土丘陵区不同土地利用下土壤粒径分析中表明, 多重分形维数与有机碳显著相关, 可以为描述土壤基本性质提供科学参考。因此, 利用分形模型研究土壤的粒径分布对于理解当地草地恢复的演变过程, 以及评价当地土壤质量有着至关重要的作用。前人大都使用单重分形描述土壤粒径分布^[14], 但其只能描述土壤结构的整体性特征^[15], 不能表征具体粒径的分布范围。而多重分形维数能详细反映土壤粒径分布的均匀程度, 并且表征土壤粒径分布的异质性和自相似性^[16], 可以更准确的分析土壤结构特性。因此, 二者的联合使用更能全面揭示土壤粒径的分布特征, 但目前的研究多以单一分形或多重分形为主, 将两种方法结合的研究较少, 尤其是生态环境比较脆弱的青藏高原。鉴于此, 本研究在青藏高原选取了4个不同利用方式的草地(放牧草地、围栏禁牧草地、围栏禁牧+补植草地和未干扰草地), 结合单重分形和多重分形理论, 研究不同草地利用对土壤粒径分布的影响, 探讨土壤粒径分布与土壤结构、土壤理化性质、土壤质量的关系, 以期为退化草地的恢复提供理论依据。

1 研究区概况

试验区(91°40′E, 32°21′N)位于西藏自治区林芝市工布达江县邦杰塘草原野外试验场。该区属于高原温带半湿润气候, 地势平坦, 平均海拔4672 m, 年平均气温-3.1℃, 最冷1月平均气温-15℃, 最热月7月平均气温9.3℃。年平均降雨量409 mm, 主要集中在夏季, 风速≥17 m/s的天数为139 d。年积温846℃, 年均日照2580 h, 无绝对无霜期。该区土壤类型为高山草甸土, 土壤质地为砂质土(图 1), 植被类型为典型的高寒草地, 主要植被有高山蒿草(Kobresia pygmaea)、紫花针茅(Kobresia pygmaea)、马先蒿(Pedicularis reaupinanta)、高山委陵菜(Potentilla polyschista)、垂穗披碱草(Elymus nutans)。

2 材料与方法 2.1 实验设计与土壤样品采集

采样区面积约为200 hm², 在研究区内设置了4种草地管理类型: 放牧草地、围栏禁牧草地、围栏禁牧+补植草地以及未干扰草地。试验区地势平坦, 无水蚀现象发生, 且地表无风蚀造成上凸下凹及地表粗化现象。所选各处理植被特征如表 1所示。每种处理5个重复小区, 共20个试验小区, 每个小区占地1 hm²。处理方法如下: (1)放牧草地: 每公顷每年放牧4头牦牛, 其放牧头数根据研究区的平均生物量、牦牛进食量、放牧面积以及放牧时间来确定。(2)围栏禁牧草地: 2009年, 采用围栏禁牧方式防止食草动物进入, 在围封之前该草地与放牧草地相同。(3)围栏禁牧+补植草地: 2009年, 在禁牧草地的基础上, 种植乡土优势植物披碱草(Elymus nutans griseb)、冷地早熟禾(Poa cymphila)和梭罗草(Kengyilia thoroldiana), 播种密度分别为3、5和8 g/m²的速率播种, 播种穴深3 cm, 间距20 cm。(4)未干扰草地: 1999年之前为放牧地, 1999年之后完全禁牧, 至今无人为干扰以及放牧活动的草地。

2018年9月采集土壤样品, 在各处理中的小区中沿对角线建立9个10 m×10 m的样方, 相邻样方间隔为10-20 m, 在每个样地内中随机设置5个1 m×1 m的小样方, 在样方内按“S”型采用土钻取0-15 cm土壤样品, 混合。将9个样地的土壤混合形成一个土壤样品, 即每个小区1个混合土样, 将所有土壤带回实验室风干过筛。

2.2 样品处理与分析

土壤颗粒组成使用MS2000型激光粒度分析仪进行测定, 仪器测量范围为0.02-2000 μm, 搅拌速度2500 r/min, 遮光范围10%-20%。称取0.5 g过2 mm筛的风干土壤, 加入10 mL10%的H₂O₂沙浴加热以去除土壤中的有机碳, 随即加入10 mL 10%的HCl煮沸使其反应完全。在上清液中加入六偏磷酸钠分散剂并使用超声波清洗机振荡, 随即用MS2000粒度仪进行土壤颗粒分析。

土壤理化性质采用常规方法测定^[17]: 土壤pH采用电位法测定；土壤有机碳含量用重铬酸外加热法测定; 全氮含量用凯氏定氮法测定；全磷采用NaOH熔融-钼锑抗比色法测定；使用有机碳分析仪(Phoenix, 800TOC)测定浸提液来分析水溶性有机碳和水溶性有机氮；铵态氮和硝态氮含量采用2 mol/L KCl以1 ∶ 5的比例提取新鲜土样, 浸提液使用3-AA3连续流动自动分析仪测定；速效磷含量用Olsen方法测定。土壤理化性质见表 2。

2.3 土壤分形维数 2.3.1 单重分形维数

利用激光粒度仪获取土壤颗粒体积分布数据, 采用土壤颗粒体积分形模型来计算单重分形维数^[19], 计算公式如下:

(1)

式中, R_i为某粒径区间上下限算术平均值, R_max为最大粒径(本研究中R_max=2), V_{(r < Ri)}为土壤粒径小于R_i的土壤颗粒累计体积, V_r为土壤颗粒总体积, D为土壤颗粒体积分形维数。

2.3.2 多重分形维数

多重分形是定义在分形结构上的有无穷多个标度指数所组成的一个集合, 是通过一个谱函数来描述分形结构上不同的局域条件、或分形结构在演化过程中不同层次所导致的特殊的结构行为与特征, 是从系统的局部出发来研究其整体的特征, 并借助统计物理学的方法来讨论特征参量的概率测度的分布规律^[20]。取激光粒度仪测量范围为I=[0.02, 2000], 将其划分为对数等差递增的100个小区间I_i=[φ_i, φ_i+1], i=1, 2, …, 100, 其中lg(φ_i+1/φ_i)为一个常数, 这样划分区间即遵循了土壤粒径分布规律并且更加客观。利用多重分形方法分析区间I的土壤颗粒粒径分布特征, 须使各子区间长度相同, 令ϕ_i= lg(φ_i, φ₁), 构造一个新的无量纲区间J =[lg (0.02 /0.02), lg(2 000 /0.02)]=[0, 5], 并且生成了100个新的等距离小区间J_i=[ϕ_i, ϕ_i+1], i=1, 2, …, 100。用v表示100个区间粒径对应的体积分数, 即v₁, v₂, …, v₁₀₀。

对土壤粒径分布进行多重分形分析, 用N个尺度为ε的相等的“盒子”划分整个土壤粒径分布跨度。在区间J中, ε=5 × 2^-k, N= 2^k, 本文k取值为[1, 6], ε取值为2.5、1.25、0.625、0.312、0.156和0.078 ^[21]。μ_i (ε) 为每个子区间土壤粒径分布的概率密度, 即子区间J_i内所有的测量值V_i的加和, 其中, i=1, 2, …, 100。利用μ_iε构建一个配分函数族为

式中q为实数; μ_i(q, ε) 为第i个子区间q阶概率, 则粒径分布的多重分形广义维数谱为

(2)

(3)

利用式(2)可以得到广义维数谱D(q), 当q=0时, D₀为信息维数；q=1时, D₁为信息维数；q=2时, D₂为关联维数。

2.4 数据处理

采用Pearson相关分析和逐步回归分析评价土壤颗粒组成、土壤分形维数与土壤理化性质的关系。采用冗余分析(RDA)对分形维数变化的环境驱动因素进行分析。本研究中响应变量包含D、D₀、D₁、D₁/D₀、D₂, 解释变量为速效磷、铵态氮、硝态氮、有机碳、全氮、溶解性有机碳、溶解性有机氮、全磷、粉粒、黏粒、砂粒、pH, 包含共12个指标。所有统计分析和作图在Rv.3.6.0中进行, 图表中数据为平均值±标准误。

3 结果 3.1 不同草地利用方式下土壤颗粒分布特征

本研究中粒径分类采用美国制分类标准, 从图 2可知, 不同草地利用方式的土壤颗粒体积分数差异显著(P＜0.05), 主要以粉粒、极细砂粒和细砂粒为主。其中, 砂粒的体积分数最高(47.37%-76.95%), 其次为粉粒(16.26%-40.81%), 黏粒的体积分数最低(6.78%-12.98%)。黏粒体积分数表现为围栏禁牧草地和未干扰草地高于放牧草地和围栏禁牧+补植草地；粉粒体积分数依次为放牧草地 < 围栏禁牧+补植草地 < 围栏禁牧草地 < 未干扰草地, 分别为40.81%、29.28%、20.40%、16.26%。相较与放牧草地和围栏禁牧+补植草地, 围栏禁牧草地黏粒、粉粒含量显著增加(P＜0.05), 增幅分别为60.0%-91.4%、43.5%-80.1%。极细砂粒体积分数未干扰草地最高, 放牧草地最低。细砂粒、中砂粒、粗砂粒体积分数依次为未干扰草地 < 围栏禁牧草地 < 围栏禁牧+补植草地 < 放牧草地。

土壤颗粒分布频率曲线可以反映土壤质地情况, 曲线的变化幅度越小, 颗粒异质程度越大^[13]。图 3为4种不同草地利用方式下土壤的颗粒分布频率曲线, 未干扰草地和禁牧草地曲线变化幅度较小, 各粒径段土壤体积含量分布相对均匀, 分布异质性较大, 细颗粒(< 50 μm)含量较高；禁牧+补植草地和放牧草地的曲线变化幅度较大, 颗粒主要集中在50-500 μm粒径范围, 以粗颗粒为主, 异质性较小, 非均匀程度小。

3.2 不同草地利用方式下土壤单重分形维数

利用式1, 通过最小二乘法拟合回归曲线, 计算出4种草地利用方式下土壤颗粒的单重分形维数, 如表 3所示, 4种草地拟合方程的决定系数均接近1, 说明线性拟合效果较好。从不同草地利用方式来看, 体积分形维数依次为放牧草地 < 围栏禁牧+补植草地 < 未干扰草地=围栏禁牧草地, 其中, 围封禁牧草地单重分形维数D最高；放牧的单重分形维数最小(2.56), 质地较粗。

图 4为土壤颗粒组成与单重分形维数D回归分析结果, 分形维数D与粉粒、黏粒体积含量呈极显著的正相关关系, 与砂粒相关性较差。从关系系数R²来看, 黏粒含量与分形维数D的相关系数较大, 粉粒次之, 砂粒最小, 说明土壤分形维数和黏粒的体积含量之间相关性最大。

3.3 不同草地利用方式下土壤多重分形维数

取q值变化区间[-10, 10], 根据公式(2)、(3), 计算出4种不同草地利用方式下广义维数普D(q)。由图 5可知, 4种不同草地利用方式D(q)-q曲线呈现反“S”型单调递减, 并且具有一定的宽度。土壤分形越均匀, D(q)-q越趋近一条直线；而在非均匀分形下, D(q)变化范围较大, 土壤粒径呈现非均匀分布, 因此进一步说明了对土壤粒径进行多重分形分析的重要性。q < 0时D(q)的下降趋势比q>0显著, 说明D(q)在稀疏区域更为敏感。

表 4为不同草地利用方式下土壤广义维数谱维数。容量维数D₀可以反映土壤粒径分布范围, D₀值越大表明粒径分布范围越广, 4种草地下D₀值变化差异不显著。从信息维数D₁来看, 不同草地利用方式的D₁值介于0.86-0.89之间, 表现为放牧草地 < 围栏禁牧+补植草地 < 围栏禁牧草地 < 未退化草地。其中, 未退化的D₁值为0.89, 表明各粒径分布不规律, 异质性较大；放牧草地最低, 粒径分布异质性较小。D₁/D₀可以衡量土壤粒径分布的集中程度^[22], 4种草地依次为放牧 < 围栏禁牧+补植草地 < 围栏禁牧草地 < 未退化草地, 未退化草地土壤粒径分布最集中, 放牧草地次之, 放牧草地最分散。关联维数D₂反映土壤粒径的均匀性, D₂值越大, 分布越均匀, 未退化草地、围栏禁牧+补植草地、围栏禁牧草地、放牧草地的值分别为0.88、0.82、0.85和0.82, 其中, 未退化草地和放牧草地的D₂值最高, 土壤颗粒各分级百分比趋向一致；放牧最低, 土壤颗粒分布均匀性最低。

3.4 土壤理化性质、土壤颗粒组成以及分形维数之间的关系

不同草地利用下土壤理化性质、土壤质地以及分形维数的相关分析如表 5所示, 黏粒与pH呈显著正相关(P＜0.05)。粉粒与信息维数D₁与有机碳、全氮和硝态氮表现极显著正相关(P＜0.01), 与全磷表现显著相关(P＜0.05), 其中与有机碳相关性最大, 相关系数为0.647。砂粒与土壤理化性质表现负相关关系, 其中与全氮呈现极显著负相关(P＜0.01), 相关系数为-0.707。单重分形维数D与土壤理化性质没有相关性, 多重分形维数(D₁、D₁/D₀、D₂)与有机碳、全氮、全磷和硝态氮呈显著正相关(P＜0.05)。

从图 6可知, 土壤基本性质共解释了96.78%分形维数的变异, 其中第一轴的解释率为71.44%, 第二轴的解释率为25.34%。砂粒、黏粒、SOC、TN和粉粒对D、D₀和D₂解释率较大。砂粒对D₀的影响最大, SOC对D₂的影响最大, 铵态氮对D的影响最大。逐步回归结果显示(表 6), 粉粒和黏粒能共同解释D值99%的变异, 表明黏粒和粉粒是不同草地利用方式下D值变化的主控因素；粉粒能够独立解释D₁59.85%的变异, 是D₁变化的主控因素；砂粒能够分别解释D₁/D₀、D₂值82.2%和79.5%的变异, 表明黏粒和粉粒是不同草地利用方式下D₁/D₀、D₂值变化的主控因素。

4 讨论 4.1 不同草地利用方式对土壤颗粒分布特征的影响

土壤颗粒组成受母质特征以及环境变化的影响^[23], 在一定程度上可以决定土壤的基本性状。研究区位于青藏高原高寒气候带, 由于气候条件影响, 草地生态系统结构单一, 在环境以及人为活动干扰下草地退化日趋严重^[24]。草原气候干燥, 在风力侵蚀作用下, 细颗粒流失, 土壤逐渐粗粒化。在本研究中, 4种草地土壤颗粒主要分布在20-250 μm之间, 0-2 μm和500-2000 μm范围内粒径含量较少, 林永崇^[25]在对藏北高原土壤的研究中也发现了类似的结果, 这是由于样地位于荒漠化地带, 成土母质源于风尘沉积, 土质疏松, 以风沙土为主。有研究表明, 围封禁牧可以有效促进土壤颗粒细化^[26]。本研究发现, 放牧草地和围栏禁牧+补植草地的砂粒含量较多, 细颗粒较少；围栏禁牧草地和未干扰草地的黏粒与粉粒含量相对较多。谢莉等^[27]研究发现, 相较于放牧草地, 围栏封育草地中土壤细颗粒含量显著增加。这可能是由于采取禁牧措施后, 植被生长恢复, 覆盖度提高, 风蚀作用减弱, 减少了细颗粒的流失, 同时植被能够截存降尘^[26], 外来的风蚀物质依附在植物表面随之沉积, 成为土壤细颗粒重要组成部分。从土壤颗粒分布频率曲线(图 2)来看, 未干扰草地变化幅度最小, 质地较均匀, 这说明未干扰草地由于植被覆盖以及矿物质含量的供应, 土壤发育良好, 且地下根系较多, 须根发达, 能够固结土体、改善土壤结构。

4.2 不同草地利用方式对分形维数的影响

单重分形维数可以反映土壤颗粒物质的损失状况, 是表征土壤物理性质的重要工具^[28]。本研究发现, 和未干扰草地相比, 其他草地土壤质地都逐渐粗颗粒化, 各草地间土壤粒径分布差异明显, 分形维数也有所不同, 说明不同草地利用方式是影响土壤质地的重要因素, 这与罗楠^[29]的研究相似。文海燕等^[30]研究表明, 围封禁牧后表层土壤单重分形维数显著增加；李国旗^[31]的研究结果显示围栏封育内的单重分形维数大于围栏外的单重分形维数。本研究也得到了相似的结论, 4种草地利用方式下分形维数D的分布范围为2.56-2.67, 依次为放牧草地 < 围栏禁牧+补植草地 < 未干扰草地 < 围栏禁牧草地。放牧草地和围栏禁牧+补植草地分形维数明显低于围栏禁牧草地和未干扰草地。未干扰草地植被常年覆盖, 土壤发育能力强, 其土壤粒度组成越细；实行围封禁牧后的草地, 人为活动少、干扰小, 植被覆盖度逐渐提高, 降低了风蚀的影响, 促进细粒物质沉积；经过补植后的草地虽然放牧活动减少, 但是草原环境恶劣, 植物生长速度较慢, 由于弥补损耗使得土壤结构不稳定。

本研究中4种草地利用方式下的土壤都具有多重分形特征, D(q)-q曲线呈现反“S”型单调递减, q < 0时D(q)的下降趋势比q>0显著, 说明D(q)在稀疏区域更为敏感, 这与王德^[13]研究结果相似。其中, q < 0时, 放牧草地与其他草地利用方式相比其变化趋势最小, 波普最窄, 说明放牧草地的土壤分形结构的复杂程度较低；q>0时, 未干扰草地变化趋势较小, 说明其土壤质地较均匀, 这与本文单重分形结果一致。D₀可以作为土壤质地分类的维数^[21], 本研究表明, 4种草地利用方式下D₀变化不明显, 孙梅^[12]研究中也得到了类似的结果表明4种草地利用方式并没有改变土壤质地组成。本研究表明, 4种草地利用方式下D₁、D₁/D₀和D₂均遵循放牧草地 < 围栏禁牧+补植草地 < 围栏禁牧草地 < 未干扰草地这一规律, 与围栏禁牧+补植草地和放牧草地相比, 围栏禁牧显著提高土壤粒径集中程度, 呈现不规律分布, 异质性较大, 土壤逐渐细颗粒化。原因可能是围栏禁牧后植被恢复、凋落物增加, 提高了土壤有机质的含量, 有机质利于土壤细颗粒的形成；青藏高原退化过程较为复杂, 除了人为活动的干扰还受自然环境的制约, 同时由于实验处理的年限限制, 围栏禁牧+补植草地植物凋落物长时间堆积, 影响周转效率, 不利于草地土壤稳定。

4.3 分形维数与土壤颗粒组成与理化性质的关系

土壤颗粒组成能够影响土壤养分的供应能力, 对植被的生长也有一定的促进作用^[32]。本研究中, 有机碳、全氮和全磷与粉粒呈显著正相关, 罗雅曦^[33]研究中也得到了相似的结果, 一方面随着植被演替的进行, 地上生物量逐渐增多, 凋落物在微生物的作用下不断分解, 导致土壤碳、氮含量增加, 另一方面, 植被能有效减少土壤细颗粒的吹蚀, 同时细颗粒具有粘结性, 使得土壤中的碳、氮含量得到了保持。在荒漠化地带, 黏粉粒的粘结作用能够促进土壤-植被系统的良性循环, 这也说明土壤颗粒组成对有机碳的固持有至关重要的作用^[34]。

土壤颗粒组成是影响分形维数的重要因素之一。已有研究表明, 土壤颗粒组成与单重分形维数有显著的相关关系^[19-35], 本研究发现, 单重分形维数D与黏粒体积含量呈极显著正相关, 与粉粒呈显著正相关(图 4)。魏茂宏^[36]在高寒草甸退化土壤粒径分布的研究中发现, 单重分形维数与黏粒百分含量呈显著正相关, 说明单重分形受土壤细颗粒含量的影响, 能够表示土壤颗粒组成情况并且反映其变化特征, 土壤颗粒越小, 分形维数越高。以上表明利用单重分形维数描述高寒高原土壤的质地变化具有一定的可行性, 可以用其来表征荒漠土壤随环境变化的趋势, 并且能够作为土壤质量演变的评价指标。而王瑞东^[37]通过对希拉穆仁草地的土壤分形维数进行研究发现, 该地区单重分形维数的大小主要与土壤粉粒含量有关, 黏粒和砂粒含量对其影响较小。可能是由于研究区的环境和土壤的差异, 同时, 在采样过程中, 取样的随机性也会对单重分形维数产生影响。

张毅^[38]在研究黄河三角洲地区土壤颗粒的分形特征中发现, 单重分形维数D与有机质呈显著正相关；杜雅仙^[39]的研究表明单重分形维数与有机质和全氮含量呈显著负相关, 与全磷和全钾含量关系不显著, 而本研究发现, 分形维数D与土壤理化性质显著关系不明显, 这与单桂梅^[40]的研究结果一致, 表明土壤单重分形在表征土壤质量上具有加强的空间特异性。相较于D值, D₁、D₁/D₀和D₂与土壤理化性质具有较好的相关性, 因此多重分形维数能更好的反映土壤质量特征^[41]。逐步回归分析表明(表 6), 粉粒是不同草地利用方式下D、D₁值变化的主控因素, 这与白一茹^[19]的研究相似, 说明随着粉粒含量增多, 土壤颗粒分布呈现异质性和不规律性。而孙梅^[12]在研究红壤多重分形特征中发现, 多重分形维数与粉粒呈负相关关系, 这可能与研究区的土壤质地有关, 该地为红壤土, 而本研究土壤以砂质土为主, 由于研究区质地不同, 导致多重分形维数也因此不同。在本研究中, 多重分形维数(D₁、D₁/D₀与D₂)和土壤理化性质表现极显著正相关, 孙哲^[42]在高寒草甸退化的研究发现, 多重分形维数与有机质含量呈显著相关。土壤有机质是土壤质量评价的重要指标^[43], 有机质含量与D₁, D₂, 和D₁ /D₀的相关关系使得这些分形维数可能作为反映土壤质量的潜在指标, 说明多重分形维数可以为衡量土壤结构和质量提供参考。

5 结论

单重和多重分形结果表明, 禁牧促使土壤粒径分布异质性更大, 颗粒分布更加集中, 质地更均匀, 且提高了细颗粒含量, 促进砂粒向黏粒和粉粒转变, 改善了土壤质地和结构, 可以作为草地恢复的参考, 但是与未干扰草地相比, 还有一定差距。单重分形维数与黏粒、粉粒呈极显著正相关(P＜0.01), 与沙粒呈负相关(P＜0.01), 可以用来反映土壤的物理性质, 而多重分形维数与土壤有机碳、全磷、全氮表现极显著的正相关关系, 可以作为衡量土壤养分的重要指标。