生态化学计量学是分析生态系统能量、元素(主要为C、N和P)组成平衡及耦合关系的综合型科学，是研究多种化学元素的平衡状态对生态交互作用影响的新理论^[1]。最初，生态化学计量学的研究主要面向植物、植被组织^[2]，后来延伸到土壤方面，并在土壤养分循环及限制作用的研究与应用中受到了较多关注^[3]，也为研究干旱瘠薄区不同生态修复措施实施效果及限制因素提供了新方向。研究区位于黄河一级支流10大孔兑之一的罕台川孔兑合同沟小流域，地处达拉特旗境内的砒砂岩裸露区；合同沟小流域是黄河上中游管理局设定的重点治理小流域，属于典型的砒砂岩强烈侵蚀区，流域内水土流失剧烈，是黄河多沙粗沙的主要策源地之一，也是达拉特旗乃至鄂尔多斯多沙粗沙的代表区域^[4]。为治理区域内剧烈的水土流失、修复生态环境，该区于20世纪90年代被列为世行贷款第1期小流域后，流域的治理趋于规范化和科学化，此后小流域内实施了一系列生态退耕造林还草工程，使当地生态环境及水土流失情况得到了一定改善，但砒砂岩地区整体生态环境依旧比较脆弱，水土流失、荒漠化问题依然存在^[5]；因此，通过研究不同退耕还林措施的成效，可为选择切实有效的生态修复措施提供依据，对治理该区的水土流失、恢复流域植被和建设生态环境有重要的现实意义。

土壤颗粒的大小能够显著影响土壤有机质的亲和性、离子的交换和植物生理活动^[6]，生态修复措施可以增加土壤黏、粉粒含量的比例，改变土壤粒径分布特征^[7]，进而改善土壤质地及养分含量，因而土壤粒径分布(particle size distribution，PSD)的变化会对土壤生态化学计量特征产生重要影响。目前，分形理论是研究具有分形特性的复杂不规则几何体最有效的方法^[8]，以往研究中诸多学者使用单重分形(D_V)作为评价土壤PSD状况的依据，但D_V虽可描述土壤颗粒分布的整体性与均一性程度，却难以描述出PSD特征的复杂空间变异性及非均质行为^[9]。多重分形在土壤学中的运用弥补了这一不足，结合多重分形维数有利于探究土壤颗粒的局部空间异质性对养分元素及化学计量的影响。目前有关砒砂岩区的研究多是以土壤理化性质改良评价为主^[10]，结合土壤颗粒分布运用多重分形探讨土壤颗粒与生态化学计量特征关系的研究还尚属空白。

基于此，本研究以合同沟小流域本氏针茅(Stipa bungeana Trin.)草地和油松(Pinusta buliformis Carr.)、柠条(Caragana korshinskii Kom.)、小叶杨(Populus simonii Carr.)退耕还林地为研究对象，以荞麦(Fagopyrum esculentum Moench.)坡耕地为对照，探究不同退耕还林措施下PSD及土壤养分含量分布的差异性及其原因，为砒砂岩区退耕还林工程改良成效、限制因素评价及树种的合理选择提供理论依据。

1 研究区概况

研究区地处内蒙古达拉特旗，属于黄河一级支流10大孔兑之一的罕台川孔兑合同沟小流域，地理位置为E 109°53′36″~110°06′53″，N 39°59′58″~40°13′18″。属温带半干旱大陆性气候，多年平均降水量310.3 mm，71.2%的降水集中在7—9月份，平均蒸发量2 157 mm。研究区为典型的砒砂岩裸露区，是达拉特旗乃至鄂尔多斯多沙粗沙的代表区域，区域内岩性结构松散，抗蚀能力差，水土流失严重，土壤质地为砂质。区内主要建群植物有羊草(Leymus chinensis(Trin.) Tzvel.)、本氏针茅等。水土保持树种主要有柠条、油松、小叶杨、沙棘(Hippophae rhamnoides L.)等。

2 材料与方法 2.1 样地选取

研究区自1998年起在区域内实施了大面积退耕还林措施，将原耕地整地后种植油松、柠条、小叶杨1~2年生幼苗，部分坡耕地弃耕后使其自然演变，并同步采取封育禁牧措施，减少人为扰动。2018年7—8月对小流域进行实地调查后发现，1998年退耕还林工程项目区平均坡度介于5°~15°之间，造林主要树种油松、柠条、小叶杨幼苗均已成林，弃耕地演变为该区地带性本氏针茅草地，而该区现有坡耕地主要以荞麦地为主^[11]。在查阅小流域实施退耕还林的设计中得知，退耕还林前土壤质地为砂质。为使采样地点在研究区有代表性，不同退耕还林的林地间有可比性，在全面踏勘的基础上，结合小流域内退耕还林工程实施时营造的水土保持林及封育的弃耕地，本研究选择坡面较平缓、坡度相近、坡向为半阴坡、光照条件相近且树种长势良好、植被盖度约30%~34%、1998年的项目区内不同树种的林地为样地类型，不同样地类型主要为以油松(YS)、柠条(NT)、小叶杨(XY)等为主的水土保持林地和本氏针茅草地(ZM)，并以当地仍然耕作的荞麦坡耕地(QM)为对照。各调查样地林木株数均>300株，各调查样地概况见表 1，各样地位置见图 1。

在各样地类型的坡面上，沿对角线分别设置3块地势相对平缓、坡向为半阴坡、坡度、光照条件植物长势基本一致的不同树种林地坡面上，按照生态学基本调查方法，设置20 m×20 m的样地，在各样地内以5点取样法取样，在各样地开挖0~60 cm土壤剖面，以10 cm为1层，共取6层，采集样品数450个。取样后，将样品带回实验室，将同一样地内的同一土层深度土壤混合，混合后样品数为90个，去除植物根系及残体，风干后研磨过2.000和0.149 mm筛，用于土壤粒径及土壤C、N和P含量的测定。

2.2 样品测定与数据分析 2.2.1 样品测定

土壤粒径采用马尔文激光粒度仪测定；土壤全氮采用凯氏定氮仪测定；有机质采用重铬酸钾-外加热法测定，经转换系数求得有机碳含量；土壤全磷采用碱融-钼锑抗比色法测定；土壤C、N和P含量取平均值，平均值与标注误差均精确至1‰。具体操作方法参照文献[12]。

2.2.2 土壤PSD多重分形参数计算

激光粒度仪测量粒度区间I=[0.02, 2 000]，以对数差划分为区间I_i=[Φ_i，Φ_i+1]，i为1~100的正整数，将区间对数变化φ_i=lg(Φ_i/Φ₁)，得到纲量一区间T=[0, 5]，内有100个等距子区间T_i=[φ_i，φ_i+1]，区间T含有N(δ)=2^k个相同尺度的小区间，其尺寸为δ=5×2^-k，k取值范围为1~6，P_i(δ)为各子区间土壤粒径分布概率密度(%)，利用尺度δ、参数q和P_i(δ)构造分配函数族^[9]：

$ {u_i}\left( {q, \delta } \right) = \frac{{{P_i}{{\left( \delta \right)}^q}}}{{\sum\limits_{i = 1}^{N\left( \delta \right)} {{P_i}{{\left( \delta \right)}^q}} }} \circ $

(1)

式中：u_i(q，δ)为第i个子区间的q阶概率；q为实数。

PSD的多重分形广义维数D(q)：

$ D\left( q \right) = \frac{1}{{q - 1}}\mathop {\lim }\limits_{\delta \to 0} \frac{{\log \sum\limits_{{\rm{ }}i = 1}^n {P_i^q} \left( \delta \right)}}{{\log \left( \delta \right)}}\left( {q \ne 1} \right), $

(2)

$ D\left( q \right) = \mathop {\lim }\limits_{\delta \to 0} \frac{{\sum\limits_{i - 1}^n {{P_i}\left( \delta \right)\log {P_i}\left( \delta \right)} }}{{log\left( \delta \right)}}{\left( {q = 1} \right)_ \circ } $

(3)

PSD的多重分形奇异谱指数：

$ \alpha \left( q \right) = \mathop {\lim }\limits_{\delta \to 0} \frac{{\sum\limits_{i = 1}^n {{\mu _i}\left( {q, \delta } \right)\log {P_i}\left( \delta \right)} }}{{\log \left( \delta \right)}} \circ $

(4)

相对于α(q)的PSD多重分形谱函数：

$ f\left( q \right) = \mathop {\lim }\limits_{\delta \to 0} \frac{{\sum\limits_{i = 1}^n {{\mu _i}\left( {q, \delta } \right)\log {\mu _i}\left( {q, \delta } \right)} }}{{\log \left( \delta \right)}} \circ $

(5)

q取值范围为[-10, 10]，步长为1，通过最小二乘法对式(2)~(5)拟合得D(q)、α(q)和f(q)。关于D(q)，本文只讨论q为0、1时的特殊情况，即容量维数(D₀)和信息维数(D₁)。

2.3 数据处理

试验数据采用EXCEL 2016进行计算处理，采用SPSS 24.0软件对各指标进行统计分析。采用最小显著差异法(least significance difference, LSD)进行差异显著性比较，显著水平为显著(P < 0.05)和极显著(P < 0.01)。采用Origin 2017进行图件的绘制。

3 结果分析 3.1 不同退耕还林措施下的PSD特征

不同措施下PSD特征有所差异(表 2)，退耕还林后土壤的细、中砂粒(0.10~0.50 mm)细化为黏、粉粒(0~0.05 mm)，柠条、油松、小叶杨和本氏针茅草地的土壤粉粒含量分别比荞麦坡耕地增加83.59%、30.14%、10.00%和29.43%，土壤砂粒含量分别比坡耕地减少36.3%、12.9%、4.5%和12.6%。除柠条林以外，其余措施对黏粒的增加效果不明显。

谱宽Δα(α_max-α_min)和谱形Δf(f(α)_min-f(α)_max)是刻画土壤分型特征的重要参数，其中：Δα表示土壤分形测度的均匀度，Δα值越大表明土壤PSD分布的不均匀性越高；Δf表示不同粒度级所处的地位，Δf＞0时，图像呈左钩上凸曲线，此时土壤粒径分布的大概率子集占主要地位，反之相反^[7]。不同退耕还林措施PSD分布的奇异谱指数α(q)-多重分形谱函数f(q)图像具有差异性(图 2)，Δα均值介于0.192~0.452，荞麦坡耕地最小，柠条、油松、小叶杨林地和本氏针茅草地的Δα值分别比坡耕地增大135.24%、124.35%、97.45%和59.38%。不同措施下Δf均为正数，图像呈左钩，说明研究区粒径分布的大概率子集是决定PSD特征的主要因素，Δf均值介于1.139~1.721之间，表现为柠条林＞油松林＞小叶杨林＞本氏针茅草地＞荞麦坡耕地。由上述结果可知，退耕还林后不同措施PSD分布的不均匀性均有所提升，并且林地的改良效果优于草地。

D₀表征土壤PSD颗粒群体的范围，数值越大则缺失粒径越少^[9]。不同措施D₀均值介于0.848~ 0.904，大小关系为柠条林(0.904)＞油松林(0.884)＞小叶杨林(0.875)＞本氏针茅草地(0.856)＞荞麦坡耕地(0.848)(表 3)。退耕还林后土壤PSD分布范围均有所增加，其中柠条林的效果最明显，即砒砂岩区柠条灌丛增加PSD范围的效果最好(P < 0.05)。D₁表征PSD均匀度，数值越大则PSD不均匀度越高^[9]。D₁均值介于0.753~0.774之间，荞麦坡耕地D₁值最小，分别是柠条、油松、小叶杨林和本氏针茅草地的97.3%、98.8%、99.3%和99.0%。由此，退耕还林措施可以提高PSD的非均一性，其中柠条林的改良效果显著(P < 0.05)，其余措施的改良效果不明显。D₁/D₀表征PSD的异质程度，数值越大则离散程度越小^[9]。D₁/D₀均值范围为0.856~0.892，不同措施及坡耕地大小关系与D₀均值完全相反，大小关系为荞麦坡耕地(0.892)＞本氏针茅草地(0.888)＞小叶杨林(0.866)＞油松林(0.862)＞柠条林(0.856)，说明退耕还林有利于增大PSD的异质性，使得不同粒径的颗粒在土壤中均有分布，林地对土壤PSD离散程度的改良效果优于草地，其中柠条林改良效果最为显著(P < 0.05)。

3.2 不同退耕还林措施下土壤C、N、P元素含量及其计量特征 3.2.1 土壤C、N、P含量

实施退耕还林后土壤C、N、P元素含量均显著高于荞麦坡耕地，但整体含量依然较低(图 3)。土壤SOC(soil organic carbon)平均质量分数为3.323~5.948 g/kg，柠条、油松、小叶杨林地和本氏针茅草地分别比荞麦坡耕地增加了79.00%、63.49%、25.19%和26.80%。草地SOC含量随剖面深度增加而减小，而人工林呈现随剖面深度增加SOC含量先增大后减小。

土壤TN(total nitrogen)含量的变化趋势与SOC含量相同，不同措施土壤TN平均质量分数为0.369~0.696 g/kg，柠条、油松、小叶杨林地和草地土壤TN含量分别比坡耕地增加了88.62%、69.30%、40.59%和15.26%。与C素的剖面分布相同，草地表层土壤(0~20 cm)TN含量较高，人工林深层土壤(>20~40 cm)TN含量较高。不同措施间土壤TP(total phosphorus)平均质量分数的变化较小(0.040~0.063 g/kg)，大小关系为柠条林＞油松林＞小叶杨＞本氏针茅草地＞荞麦坡耕地。不同措施间土壤TP含量与SOC、TN含量呈现出较一致的剖面分布规律。

3.2.2 土壤C、N、P计量特征

不同措施间土壤C、N、P化学计量比总体变化幅度较小(图 4)。土壤C/N均值介于8.02~9.91，不同措施间C/N(carbon-nitrogen ratio)波动较小，且随剖面深度增加C/N差异不显著。土壤C/P(carbon-phosphorus ratio)均值为77.65~94.34，N/P(nitrogen-phosphorus)均值为8.50~11.03。土壤C/P和N/P随剖面深度的增加先增大后减少，林地深层土壤(30~60 cm)的C/P和N/P显著(P＜0.05)高于表层土壤(0~30 cm)，其中柠条、小叶杨林地增幅最大，深层土壤C/P和N/P分别是表层土壤的1.51、1.43倍和1.50、1.52倍。

3.3 土壤PSD分布和C、N、P元素及其计量比的相关关系

土壤PSD分布与C、N、P元素含量及其计量间的Pearson相关分析如表 4所示。SOC、TN含量与黏粒、粉粒含量呈极显著或显著正相关(P < 0.01或P < 0.05)，土壤P含量仅与黏粒含量呈显著正相关(P < 0.05)；细砂含量与SOC、TN、TP含量均呈极显著负相关(P < 0.01)。土壤SOC、TN含量与Δα、D₀呈极显著正相关关系(P < 0.01)，TP含量与分形参数间未达到显著水平。土壤C、N、P化学计量比未与土壤黏、粉、砂粒含量呈现统一的相关关系，且相互间关系不显著，C/P、N/P与Δα、Δf均呈极显著正相关(P < 0.01)，C/P与D₀呈显著性相关(P < 0.05)，其余指标间的相关关系未达到显著水平。

由此可知：黏、粉粒和细砂是决定本区域土壤C、N、P元素含量的主要土壤粒径，PSD组成及分形参数可以在一定程度上表征土壤C、N、P元素的含量；而PSD分布对C、N、P化学计量比的影响较小，谱宽(Δα)和谱形(Δf)可以较好的解释土壤的C/P和N/P特征。

4 讨论 4.1 不同措施对土壤PSD分布的影响

土壤PSD分布与土壤质地、养分持有量有着密切的联系^[8]，以往研究发现，退耕还林措施可以显著改变土壤PSD分布特征，细化土壤颗粒改善土壤结构，提高土壤D₀、D₁等分型参数值^[7]。本研究中不同措施间土壤容量维数D₀、谱宽Δα和谱形Δf差异显著(P＜0.05)，且均大于荞麦坡耕地，说明退耕还林措施有利于增加PSD分布的范围及非均匀性，优化土壤质地。而不同措施表层土壤(0~10 cm)谱宽(Δα)差异较小且显著低于其他土层(P＜0.05)，介于0.163~0.181之间。这可能是研究区低植被覆盖度及复杂的风水复合侵蚀环境共同作用的结果。植被覆盖可以有效降低近地表风速，保护土壤细粒物质。植被盖度作为覆被条件的重要指标，与风蚀强度有密切联系^[13]，研究区植被盖度介于30%~34%之间，较低的植被盖度造成植被对表土细粒物质的保护作用不足，使得土壤仍有继续流失的可能，造成不同措施土壤PSD分布均存在一定程度的细粒组分缺失的现象，因而PSD的非均匀程度被降低。同时，研究区复杂的风水蚀复合侵蚀环境也在一定程度上促进了表土细粒物质的流失，脱登峰等^[14]研究表明，水蚀与风蚀之间具有明显的正交互效应。风蚀会增强降水对坡面的侵蚀强度，而水蚀通过破坏土壤结构为风蚀提供物质来源^[15]，促进风蚀的发生。以上原因造成了表层土壤Δα值显著低于其他土层；因而本研究建议在后续恢复中采取封禁措施、人工种草措施，减少人为破坏，增加地表植被覆盖度，减少表土及细粒土壤的流失。

4.2 不同措施对土壤化学计量特征的影响

土壤C、N、P含量与植被光合作用及土壤矿质代谢等过程密切相关^[16]，并且会直接影响土壤微生物活性、枯落物分解和土壤的肥力条件^[17]，其化学计量比可以反映土壤有机构成及养分能力^[12]。本研究发现，不同措施下土壤C、N、P平均质量分数为(4.939±0.301)、(0.566±0.037)、(0.056±0.004)g/kg，远低于全国平均水平11.12、1.06、0.65 g/kg，且P含量严重匮乏，只有全国平均水平的8.65%，这一现象与本地区强烈的侵蚀环境密切相关。研究区砒砂岩土壤结构差，长期剧烈的土壤侵蚀造成土壤养分大量流失，导致C、N、P含量极低。不同人工林地深层土壤C、N含量均高于表土层，这与杜满义等^[18]研究结果有所差异。尽管土壤C、N元素的积累主要源于有机质和枯落物的分解^[19]，但研究区的风蚀、水蚀交互作用使枯落物难以在林内长期聚集，从而无法达到有机质释放SOC、TN的腐殖化阈值；深层土壤微生物对根系残体分解能够释放出大量SOC和TN^[20]，植被根系的分泌物对土壤矿质的分解作用^[21]，在一定程度上促进了根系分布密集层土壤C、N含量的增加。这也说明，在研究区退耕还林年限到达20年时，土壤养分含量均有所提高，但总体水平依然较低，表层土壤的改良效果低于深层土壤。为更好促进土壤改良效果，特别是表层土壤养分的提高，一方面可通过保护地表枯落物、增加枯落物覆盖及其水源涵养、土壤保持功能，以此减少地表土壤流失，同时利用枯落物的有效积累与分解促进养分的增加，另一方面，通过封禁、人工种草等方式促进林下植被恢复，提高植被盖度、保护表层土壤。这些措施不仅有利于促进林木的生长，而且有利于提高本地区生态修复效果。

土壤化学计量是土壤养分循环与平衡的重要参数，而C/N是影响土壤C、N循环，表征土壤质量的敏感指标^[22]。本研究发现不同措施及剖面层次间土壤C/N波动均较小。对此，杜满义等^[18]认为土壤中C、N元素间显著正相关关系以及耦合作用会降低土壤C/N变异程度，这与本研究相关分析结果相符。土壤中C、N含量间具有一定的动态平衡，当C含量增加时，土壤微生物会倾向于进行N积累过程，相对减少N素的矿化^[23]，因此C、N含量的变化会在一定范围内保持高度的一致性。

4.3 土壤颗粒对土壤化学计量特征的影响

土壤颗粒是黏结土壤有机质等营养元素的重要介质^[24]，相关分析结果表明，土壤C、N元素与黏、粉粒含量、D₀和Δα呈显著或极显著正相关(P < 0.05或P < 0.01)，说明土壤中优质细颗粒的增加有利于促进土壤C、N元素的积累，PSD分布特征在一定程度上也决定了土壤养分元素的供应能力。罗雅曦等^[6]也认为，土壤PSD分布特征及分形维数可以作为指示土壤营养元素含量的定量化指标。土壤细颗粒部分是吸附有机质的主要无机胶体，随着土壤颗粒的优化，土壤有机质合成与分解的媒介不断增加，促进了土壤C、N含量的增长。而P含量只与黏粒、细砂粒含量达到了显著水平，与其余粒径及分形参数的相关性均较低，这可能与土壤P积累主要来源于岩石的风化有关^[19]，研究区侵蚀剧烈，矿质风化程度较低，难以与植被的改良进程达到同步，因而P含量与土壤PSD分布及分形参数间相关性较低。此外，本研究发现土壤C、N、P化学计量与PSD分布及参数间未呈现统一的相关关系，且相关性均较低(Δα、Δf除外)，其原因还有待进一步研究。

5 结论

1) 砒砂岩区退耕还林措施可以将土壤的中、细砂粒细化为黏、粉粒，显著增加PSD分布范围、非均匀性和离散程度，其中对粉粒的改良效果显著(P < 0.05)，对黏粒的改良效果不明显。草地有利于表层土壤质地的细化，林地更有利于深层土壤的改良。

2) 研究区土壤C、N、P含量较低，P元素含量严重匮乏，退耕还林后土壤C、N元素含量显著提升，P素含量提升较少，不同措施土壤C、N、P含量的改良效果表现为柠条林＞油松林＞小叶杨林＞本氏针茅草地；不同措施及剖面层次间土壤C、N、P化学计量比波动均较小。

3) 黏、粉粒和细砂粒是决定本区域土壤C、N、P元素含量的主要土壤粒径，与C、N元素相比土壤P与PSD分布的相关性较弱；Δα、Δf与C/P、N/P均呈极显著正相关关系，PSD分布及参数能够较好的表征土壤C、N、P元素含量及化学计量特征。