生物结皮（Biocrusts）是由苔藓、藻类、地衣以及土壤微生物与表层土壤颗粒通过假根、菌丝体和分泌物胶结形成的复合体，覆盖了地球陆地面积的12%，在低生产率的干旱生态系统中发挥着重要作用^[1-4]。已有的研究表明生物结皮会对土壤的养分积累过程、微生物群落以及地表微地形等产生影响^[5-9]。

生物结皮不同组分的固碳、固氮作用积累的养分是旱区生态系统重要的物质输入。然而生物结皮个体矮小，且不存在真正的根系，其养分富集效应主要表现在表层土壤。Chamizo等^[5]的研究表明，生物结皮下伏土壤有机碳的含量随生物结皮的发育年限增加，Huang等^[6]发现苔藓结皮和藻结皮的固碳速率分别可达64.9和38.6 g·m^–2·a^–1，而藻结皮、地衣结皮的固氮活性可达9.81×10³ nmol m^–2·h^–1（以C₂H₄计）^[9]。目前有关生物结皮养分效应的研究多集中在地块尺度，鲜有从区域尺度上探讨生物结皮养分空间特征的报道，阻碍了生物结皮在旱区生态系统中的价值评估。

1959年以来开展的生态修复工程使毛乌素植被覆盖率达到40%~50%，流沙得到有效治理^[10]，生物结皮作为流动沙地固定的标志在该区域广泛分布，2000年以来在地块尺度的研究表明，生物结皮在该地区发挥重要的养分富集^[11-12]、水土保持作用^[13]。本文以毛乌素沙地为研究区，基于全域空间的146个样点的样品采集与指标测算，通过地统计学方法、多元统计方法等，分析生物结皮下伏土壤有机碳（SOC）、全氮（TN）、全磷（TP）含量及其空间异质性，探讨影响生物结皮下伏土壤养分变化的因素，旨在为区域尺度生物结皮养分特征的研究方法提供借鉴，为科学评价生物结皮在荒漠生态系统中的生态价值提供参考。

毛乌素沙地位于陕西省、内蒙古自治区、宁夏回族自治区三省区交汇处（37°20′~39°23′N、107°20′~111°30′E），面积为4.22万km²，沙区地表包含多种风沙地貌类型，固定沙丘、半固定沙丘和流动沙丘广泛分布，总体地势自西北向东南倾低，中西部为梁地与滩地组成的波浪高原，海拔在1 300 ~1 584 m之间，中南部地区风沙堆积增厚，向东南逐渐过渡至黄土高原的梁峁丘陵区，海拔在902 ~1 300 m之间。该区域属温带干旱、半干旱的大陆性季风气候，年平均气温为6.0℃~9.0℃，区域内水热平衡状况自西北向东南表现出有规律的变化：沙区10℃积温在2 700℃~3 000℃之间，由西北向东南逐步增加；潜在蒸散率（PER）由西北的2.5 mm·a^–1递减为东南部的0.8 mm·a^–1；太阳辐射总量自西北部的150 kcal·cm^–2·a^–1减至东南部的138 kcal·cm^–2·a^–1，年均降水量自西北向东南由200 mm递增至490 mm，7-9月降水量约占全年降水量的70%以上。区域内土壤以风沙土为主，养分含量低、抗侵蚀性能差；植被以沙生植被面积最广，主要包括沙蒿（Artemisia desertorum Spreng.）、沙柳（Salix cheilophila Schneid.）、柠条（Caragana korshinskii Kom.）等，其次为草甸植被，如寸草草甸、芨芨草草甸等。

基于2016-2017年全域的多次野外调查，布设样点时首先保证其空间代表性和分布均匀性，同时考虑生物结皮分布、样点可到达性等因素，在2016年拟定106个样点，2017年增添40个样点并完成采样，合计布设146个样点（包括72个苔藓植物主导的生物结皮样点与74个藻类主导的生物结皮样点，图 1）。记录各样点生物结皮的发育类型及其立地条件，通过手持GPS（Garmin inReach SE+）获取各样点经度、纬度及海拔信息。清除结皮层表面枯落物等杂质后，在1m×1m样方内使用微型土钻按五点采样法采取生物结皮的下伏5 cm层土壤样品，五点混匀后作为1个重复，各样点均设置3次重复采样。土壤样品装入自封袋带回室内，采用TOC分析仪（TOC-L）测定其土壤有机碳（SOC）含量，全自动间断化学分析仪（CleverChem200）测定全氮（TN）、全磷（TP）含量。

图 1 研究区及样点位置 Fig. 1 Study area and sample location

选择IBM SPSS Statistics 23对生物结皮下伏土壤的SOC、TN、TP含量进行经典统计分析。试验所得数据采用样本均值加减三倍标准差（μ±3σ）识别特异值，分别用正常的最大值和最小值代替，后续计算均采用处理后的数据。非正态分布的数据及特异值会引起空间变异分析中变异函数的比例效应，导致某些空间机构模糊^[14-15]，采用K-S（Kolmogorov-Smirnov）检验法验证数据正态性。

应用GS+ 9.0进行空间变异分析及半方差函数分析，建立最佳拟合模型后将其参数导入ArcGIS 10.3软件的Geostatistical Analyst模块中，采用普通克里格法（Ordinary Kriging）获得毛乌素沙地生物结皮下伏土壤养分的空间分布图^[16-17]。通过R软件（R-4.0.2）进行生物结皮下伏土壤的养分含量与其可能的影响因子（表 1）间的相关性分析、逐步回归分析，影响因子中的气候因子数据来源于TerraClimate数据集^[18]（1958-2015年全球月度高空间分辨率（1/24°，~4 km）气候数据集），数据获取地址为www.climatologylab.org/datasets.html；植被因子数据来自美国国家航空航天局（NASA）的地球观测系统MODIS数据集，数据通过https://ladsweb.modaps.eosdis.nasa.gov获取。

表 1 Table 1

表 1 解释变量表 Table 1 Table of explanatory variables 解释变量 注解 Explanatory variables Annotation 风速Wind speed（ws） - 降水量Precipitation（pr） - 土壤湿度Soil moisture（sm） - 平均气温Mean temperature（Tas） - 帕默尔干旱指数Palmer drought severity index（PDSI） 一种基于水量供需关系的干旱指数 饱和水汽压差Vapor pressure deficit（vpd） 表示实际空气距离水汽饱和状态的程度，即空气的干燥程度 归一化植被指数Normalized difference vegetation index（NDVI） 评价与分析植被生长、覆盖状况的主要指标之一 潜在蒸散量Potential evapotranspiration（pet） 指充分供水下垫面蒸发或蒸腾到空气中的水量，能够反映出一个地区的蒸发能力 地表向下短波辐射Downward shortwave flux at the surface（srad） 包括直射太阳辐射、大气散射辐射，是地表能量收支的重要组成部分，也是光合作用的源驱动力

表 1 解释变量表 Table 1 Table of explanatory variables

毛乌素沙地生物结皮下伏土壤的SOC含量介于0.56~20.47 g·kg^–1之间，TN含量在0.08~1.66 g·kg^–1之间，TP含量在0.17~1.18 g·kg^–1之间，研究区土壤特殊，裸沙地SOC、TN、TP的平均含量仅为0.056 g·kg^–1、0.003 g·kg^–1、0.002 g·kg^–1，生物结皮对下伏土壤三种养分含量的提升达到极显著水平（P < 0.01）（表 2）。若参照全国土壤养分含量分级标准（表 3），经换算的SOC与TN含量均处于中等水平，TP处于极低水平。变异系数（CV）直观地反映出样本数据的变异性大小，生物结皮下伏土壤养分变异程度以TP最大，其次为SOC、TN，变异系数分别为60.5%、37.9%和30.8%，均属于中等程度变异。数据分析过程中亦发现苔藓结皮下伏土壤养分含量高于藻结皮，苔藓结皮下伏土壤的SOC、TN、TP均值为11.82 g·kg^–1、0.92 g·kg^–1、0.56 g·kg^–1，藻结皮下伏土壤的SOC、TN、TP均值为9.31 g·kg^–1、0.89 g·kg^–1、0.42 g·kg^–1。

表 2 Table 2

表 2 生物结皮下伏土壤养分的描述性统计 Table 2 Descriptive statistics of soil nutrients underlying biological crusts 指标 Index /（g·kg–1） 最小值 Min. 最大值 Max. 平均值 Mean 标准误差 SD 变异系数 CV/% 偏度 Skew. 峰度 Kurt. K-S检验 K-S Test 土壤有机碳SOC 0.56 20.47 10.51 3.99 37.9 –0.02 0.17 0.08 全氮TN 0.08 1.66 0.91 0.28 30.8 –0.16 –0.26 0.20 全磷TP 0.17 1.18 0.48 0.29 60.5 0.38 –0.88 0.13（Box-Cox）

表 2 生物结皮下伏土壤养分的描述性统计 Table 2 Descriptive statistics of soil nutrients underlying biological crusts

表 3 Table 3

表 3 全国土壤养分含量分级标准表 Table 3 National soil nutrient content classification standard table 指标 Index 极高 Extremely high 高 High 中上 Upper middle 中 Medium 低 Inferior 极低 Extremely low SOM /（g·kg–1） > 40.00 30.00~40.00 20.00~30.00 10.00~20.00 6.00~10.00 < 6.00 TN /（g·kg–1） > 2.00 1.50~2.00 1.00~1.50 0.75~1.00 0.50~0.75 < 0.50 TP /（g·kg–1） > 2.00 1.50~2.00 1.00~1.50 0.75~1.00 0.50~0.75 < 0.50

表 3 全国土壤养分含量分级标准表 Table 3 National soil nutrient content classification standard table

正态性检验结果显示，SOC、TN两组数据的偏度系数分别为–0.02、–0.16，均表现出一定的左偏态；SOC数据分布相较正态分布表现为尖顶峰，TN为平顶峰，峰度系数分别为0.17、–0.26。显著性水平α=0.05的条件下应用K-S检验法所得SOC、TN的P均大于0.05，表明两组数据均服从正态分布；TP数据表现右偏态平顶分布，偏度系数、峰度系数分别为0.38、–0.88，原数据正态性较差，经Box-Cox转换后的TP含量数据服从正态分布。

地统计学分析结果（表 4）表明，生物结皮下伏土壤SOC、TN的最佳半方差函数拟合模型均为指数模型，TP为线性模型，其决定系数分别为0.816、0.742和0.652，表明各变量的拟合模型能较好地反映其空间分布特征。SOC、TP的自相关范围较大，变程分别为171.83 km和147.78 km，空间连续性较好；而TN的自相关范围较小，空间连续性相对较差。生物结皮下伏土壤SOC、TN、TP的块金系数分别为53.3%、49.86%和20.2%，SOC、TN呈中等强烈的空间自相关性，且随机因素和结构因素对其空间变异的影响非常接近，TP表现强烈的空间自相关性，其空间变异主要由结构性因素引起。

表 4 Table 4

表 4 生物结皮下伏土壤养分变异函数理论模型与相关系数 Table 4 Theoretical model and correlation coefficient of soil nutrient variation function under the biological crust 指标 Index 拟合模型 Fitting model 块金值 Nugget（C0） 基台值 Sill（C0+ C） 块基比 Nugget/ Sill C0/（C0+ C）/% 变程 Range/km 决定系数 R2 残差 RSS 土壤有机碳SOC 指数 9.171 18.35 53.30 171.83 0.816 15.6 全氮TN 指数 0.056 0.12 49.86 69.33 0.742 9.48×10–4 全磷TPBox-Cox） 线性 0.040 0.05 20.20 147.78 0.652 5.48×10–5

表 4 生物结皮下伏土壤养分变异函数理论模型与相关系数 Table 4 Theoretical model and correlation coefficient of soil nutrient variation function under the biological crust

分布规律明显、分布斑块大，是毛乌素沙地生物结皮下伏土壤养分空间分布的突出特征。三种养分的空间分布插值结果（图 2）显示，研究区生物结皮下伏土壤SOC含量在空间上自东向西呈现出高-低-高-低的带状分布特点，高值区出现在神木市西部、榆林市东北部、乌审旗西部和鄂托克旗东部，位于研究区东部及中西部，SOC含量在11.43~20.47 g·kg^–1之间，若参照全国第二次土壤普查养分分级标准，其SOC含量处于缺乏到中等含量水平；西部及中东部地区为SOC含量低值区，介于0.56~9.6 g·kg^–1之间，处于极缺乏到缺乏水平，极低值出现在横山县与榆林市之间、鄂托克前旗东部，TN与SOC含量的空间变化规律基本一致，反映在二者极显著的相关性上（r=0.818，P < 0.01，皮尔逊双尾检验），TN分布以0.47~1.08 g·kg^–1为主，处于缺乏到中等含量水平。

图 2 毛乌素沙地生物结皮下伏土壤三种养分的空间分布 Fig. 2 Spatial distribution of biological crust and its underlying soil nutrients in Mu Us Sandland

TP含量的空间分布异于SOC、TN含量的空间分布，在空间上呈现出带状与块状相结合的分布格局，西南方向出现局部的低值区，TP含量的高值区出现在研究区东部并向中部延伸，向南、北两侧其值有降低的趋势。

前文从随机因素、结构因素的角度初步揭示了研究区生物结皮下伏土壤养分的变异特征。由于研究区内地形平坦且土壤以风沙土为主，故本文关注空间异质性较强的地表水热状况、植被覆盖状况与毛乌素沙地生物结皮下伏土壤养分间的关系，并使用多元逐步回归的方法筛选主要影响因素。

通过相关系数表（表 5）可知，各气候因子以及植被因子与SOC、TN、TP三种土壤养分间的关联程度不同：3种土壤养分均与ws、srad两种气候因子表现显著的正相关关系（P < 0.05），并且SOC、TN与植被因子NDVI呈极显著正相关关系（P < 0.01）；SOC、TN、TP与Tas呈极显著负相关关系（P < 0.01），相关系数分别为–0.32、–0.28、–0.27；干旱程度直接影响生物结皮的发育状况，进而影响其下伏土壤的养分含量，PDSI亦与3种土壤养分呈负相关关系，但PDSI与TP间相关性未达到显著性水平。相关分析结果在一定程度上表明了气候因子、植被因子与生物结皮下伏土壤养分间的关系，同时也反映出解释变量间存在的多重共线性问题，比如Tas与PDSI、pet、pr、sm、srad间相关系数分别为0.55、0.59、0.33、0.45、–0.6，均达到显著性水平（P < 0.05）。

表 5 Table 5

表 5 生物结皮下伏土壤养分与影响因子间的Pearson相关性 Table 5 Pearson correlation between soil nutrients underlying biological crusts and influencing factors 环境参数 Environmental parameters ws PDSI pet pr sm srad Tas vpd NDVI SOC TN PDSI –0.58** pet –0.15** 0.34** pr –0.74** 0.42** –0.05 sm –0.79** 0.42** –0.17* 0.85** srad 0.29** –0.02 –0.04 0.17* –0.19* Tas –0.76** 0.55** 0.59** 0.33** 0.45** –0.6** vpd –0.04 0.18* 0.86** –0.25** –0.37** –0.26** 0.53** NDVI 0.02 –0.07 –0.06 0.14 0.08 0.15 –0.11 –0.09 SOC 0.24** –0.17* –0.11 0.01 –0.07 0.3** –0.32** –0.14 0.24** TN 0.21* –0.21* –0.11 –0.04 –0.1 0.23** –0.28** –0.13 0.28** 0.82** TP 0.2* –0.13 –0.02 0.08 –0.11 0.33** –0.27** –0.05 0.04 0.37** 0.28** 注：*表示P < 0.05，**表示P < 0.01，PDSI、pet、pr、sm、srad、Tas、vpd、NDVI、SOC、TN、TP分别代表帕默尔干旱指数、潜在蒸散量、降水量、土壤湿度、地表向下短波辐射、平均气温、饱和水汽压差、归一化植被指数、土壤有机碳、土壤全氮、土壤全磷。Note：* means P < 0.05，** means P < 0.01，PDSI，pet，pr，sm，srad，Tas，vpd，NDVI，SOC，TN，TP represent palmer drought severity index，potential evapotranspiration，precipitation，soil moisture，downward shortwave flux at the surface，mean temperature，vapor pressure deficit，normalized difference vegetation index，soil organic carbon，soil total nitrogen，soil total phosphorus，respectively.

表 5 生物结皮下伏土壤养分与影响因子间的Pearson相关性 Table 5 Pearson correlation between soil nutrients underlying biological crusts and influencing factors

为进一步确定气候因子及植被因子对研究区生物结皮下伏3种土壤养分变异的贡献，并处理解释变量间的多重共线性问题，采用逐步回归的方法分析SOC、TN、TP空间变异的主要影响因素，结果表明3种土壤养分的主控变异因素不同：SOC的逐步回归模型引入了土壤湿度（sm）、地表向下短波辐射（srad）、植被归一化指数（NDVI）三个变量；TN的逐步回归模型引入了平均气温（Tas）、植被归一化指数（NDVI）两个变量；风速（ws）、土壤湿度（sm）分别作为第一和第二因子被引入TP的逐步回归模型。气候因子、植被因子对研究区3种土壤养分变异的解释能力不同：对于空间自相关性强烈的TP，气候因子、植被因子能解释其空间变异的38.7%；对空间自相关性中等的SOC、TN的解释率分别为27.5%、24.1%。

与目前大多数关于生物结皮养分效应的研究结果相一致^[1-2]，本研究也发现生物结皮对下伏土壤的养分含量有显著提升。作为荒漠系统碳的主要贡献者，生物结皮所贡献的有机碳调节了土壤养分的分解和矿化率，进而调节着养分的有效性和初级生产力^[5]；Rogers和Burns^[19]发现大约70%由生物结皮固定的氮素被立即释放到土壤环境中且对周围的维管植物是有效的，是养分贫瘠的荒漠生态系统中重要的氮素来源之一。磷在土壤中的移动性弱，且生物结皮下伏土壤TP的主要影响因子同SOC、TN不同，使得TP显示出高度的变异性。SOC、TN、TP的空间异质性及其对尺度的依赖程度各不相同，表现在各养分块金系数及变程的不同：与杨之江等^[20]的研究结果类似，生物结皮下伏土壤TP表现出强烈的空间自相关性（块金系数23.2%），说明样点之间TP含量的变化主要受气候、土壤母质等因素控制；随机性因素对SOC、TN含量的影响削弱了结构性因素造成的影响，二者的空间异质性均为中等水平，这可能是因为生物结皮对C、N、P三种元素的富集能力不同，闫德仁等^[21]在对库布齐沙漠生物结皮“肥岛效应”的研究结果中发现与沙土相比，生物结皮层全氮、有机质、全磷分别提高了77.8%、55.4%、19.9%，在黄土区的研究^[22]也表明生物结皮对土壤C、N含量的影响要强于对P含量的影响。

部分研究指出，沙地生物结皮的形成在时间序列上是逐渐由物理结皮发育成生物结皮，植物冠层为下层土壤提供更多的凋落物，拦截更多的土壤微细颗粒，有利于土壤结皮的形成，另一方面维管植物提供遮荫改变地表微环境，如温度、湿度、光照强度等，为生物结皮的拓殖、演替以及养分积累过程提供了适宜的环境^[23-24]。本文的研究结果也证实了这一点：生物结皮下伏3种土壤养分的高值区在研究区东部的榆林市、神木市出现重合，SOC、TN与植被因子（NDVI）呈显著相关关系，且SOC、TN的逐步回归模型均引入了NDVI作为参数，表明这两种土壤养分的分布格局与植被覆盖状况关系密切。

除植被因子以外，土壤湿度（sm）、地表向下短波辐射（srad）分别作为第一因子和第二因子引入至SOC的逐步回归模型，这两种气候因子在研究区东部均为高值，西部（西南部）为低值，与SOC表现类似的分布状况。一方面，sm与srad通过影响生物结皮的生长发育状况进而对土壤有机碳的输入量造成影响：srad作为生物结皮光合作用的源驱动能力，是限制其生长发育的关键因子；Li等^[25]在中国北方不同降水梯度下的研究表明，生物结皮表层SOC含量与年均降水量、土壤湿度变化趋势类似，但是在本文的分析结果中，3种养分的逐步回归模型均未引入降水量（pr）作为参数，这是因为生物结皮的发育主要受近地面空气湿度影响^[26]，相较于降水量而言，土壤湿度对于生物结皮发育的影响更为显著。另一方面，适宜的水热条件有利于土壤微生物的生长繁殖^[27]，进而对SOC组分中的微生物量碳产生影响。TN的逐步回归模型在植被因子以外还引入了Tas，回归系数为–0.149，这是因为高温不利于苔藓植物生长，这一点在多数的生物结皮的培育恢复试验^[28-29]中已经得到了证实。

研究区生物结皮下伏土壤SOC、TN的空间分布相似且极显著相关。这可能与碳氮耦合作用有关，生物结皮中的苔藓植物、光合自养微生物在演替发育过程改变土壤环境条件如水分、空气、温度等，为土壤微生物提供了适宜的生长环境^[8，30]，其中的固氮微生物将大气中的氮气还原成氨，过程中消耗光合作用所提供的能量。李嵘和常瑞英^[31]的研究发现氮素积累过程会抑制参与土壤碳循环的过氧化物酶、淀粉酶等酶的活性，从而增加土壤碳积累，蔡玉婷等^[32]的研究表明氮素添加条件下土壤有酸化的风险，从而抑制土壤微生物的生长及其对土壤碳的分解。SOC、TN在研究区中东部出现条带状的低值区，这很可能受到了生物结皮类型的影响，从样点分布图（图 1）可以看出，在研究区中东部设置的采样点以藻类主导的生物结皮为主，本文的研究结果表明苔藓结皮下伏土壤养分的平均含量高于藻结皮，与目前对于不同类型生物结皮养分富集能力的研究结果相一致^[33-34]。

与SOC、TN不同，TP含量在研究区内表现出自东向西递减的趋势，其逐步回归模型分别将风速（ws）、土壤湿度（sm）作为第一因子和第二因子引入，回归系数分别为1.127、–0.059。土壤磷的扩散系数随含水量的增加而增加，而风速对研究区的大气降尘量造成影响，生物结皮（尤其是苔藓结皮）形体结构特殊，表明粗糙度大^[35]，能够捕获大气降尘中的颗粒物质使得磷素输入，闫德仁等^[36]对风沙环境下苔藓结皮层化学元素组成的研究中发现，P₂O₅在结皮层的含量明显高于对照风沙土。此外，在调查采样过程中发现，研究区西部、西南部（内蒙古自治区、榆林市靖边县）的牧场数量较多且规模大，放牧导致土壤中钾素升高，促进草本植物生长，有可能消耗生物结皮所固定的有限的磷素，研究区内生物结皮下伏土壤养分对于干扰现象的响应有待于进一步的研究。

毛乌素沙地生物结皮提升了荒漠土壤养分含量，三种养分要素均表现出明显的空间相关性，经Box-Cox变化后的TP含量表现出强烈的空间相关性。不同的拟合模型很好地反应了研究区生物结皮下伏土壤养分含量的空间特征，研究区生物结皮下伏土壤SOC、TN含量最佳拟合模型为指数模型，变程分别为171.83 km、69.33 km，决定系数分别为0.816、0.742；TP为线性模型，变程为147.8 km，决定系数为0.652。生物结皮下伏土壤养分在研究区的分布状况受到植被恢复状况与地表水热状况等因素的影响，三种养分在东部均表现为高值，其中TP含量自东向西表现递减的趋势；SOC、TN自东向西表现高-低-高-低的带状分布特点。