土壤是农业生产的基础，是人类赖以生存的基石，也是人类食物与生态环境安全的保障^[1]。土壤全碳、全氮含量可以作为评价土壤肥力的重要指标^[2]。对于土壤碳素的研究有助于了解全球碳循环趋势，进而有效地控制碳储存，减少温室气体的排放。土壤全氮是土壤氮素肥力指标，对农作物的生长具有重要作用。由于气候、植被等自然条件的多样性和人类活动等的综合影响，增加了土壤理化性质空间分异的复杂性^[3-4]，而空间异质性现象的本质为土壤空间的自相关性^[5]，研究表层土壤全碳全氮的空间自相关及环境响应机制，对深入了解大环境变化下的碳氮循环动态变化及生态修复具有不可忽视的作用。

虽有学者对黑河流域的土壤养分进行了研究^[6-8]，但针对土壤水分^[9]、有机质^[10]、盐分^[11]及呼吸速率^[12]等的研究较多，对于土壤化学因素的研究目前较少，且对于影响机制的研究多采用多元分析方法、主成分分析法及相关性分析等，这些方法仅能体现单一因素对土壤养分的影响或变量间的相关性，不能全面地反映变量之间的相关性，结果较为片面^[13]。应用地理探测器模型不仅可以全面地研究土壤养分变化的驱动因素，定量识别因子间的相互作用^[14]，而且可以研究驱动因素之间对土壤养分变化的交互作用。基于此，本研究将通过地理探测器模型探究环境因子综合变化对土壤全碳全氮含量的影响。

黑河不仅是祁连山最具代表性的河流，而且是我国西北内陆干旱区第二大内陆河^[15]。整个黑河源区不仅是祁连山国家公园的重要组成单元和重要水源补给区，而且对黑河上、中、下游地区经济社会发展和城乡居民生存环境起着决定性作用^[16]。黑河上游海拔高差大，主要以山地为主，植被类型与土壤类型丰富，土壤养分具有明显的空间异质性。在此背景下，本研究以重要的水源涵养区——黑河源区为研究区，利用地统计分析方法结合地理探测器模型对土壤全碳全氮含量进行空间变化及响应机制研究，以为该区域的生态保护及环境治理提供科学依据。

1 研究区概况与研究方法 1.1 研究区概况

研究区地处黑河上游祁连山区，地理坐标98°05′35″ ~ 101°02′06″ E，37°25′16″ ~ 39°05′18″ N，总面积近1.6×10⁴ km^2[15]。根据流向，其源头可分为东支和西支：东支源于祁连县东部峨堡镇滩东的锦阳岭，称峨堡河；西支即称为黑河，起源于铁里干山和八一冰川一带，是黑河源区的干流，途径野牛沟乡和扎麻什乡等。该区气候属于高原大陆性气候区，年降水量为340 ~ 400 mm，其中6—9月降水量占全年降水量的78% ~ 89%，年蒸发量为1 496.3 mm，年日照时数为2 500 ~ 3 100 h，太阳辐射强烈，昼夜温差大，年平均气温为–4.7 ~ 1.2℃，几乎无绝对无霜期^[16]。

1.2 土样采集与测定

综合考虑研究区的植被、海拔、土壤等自然环境因子及交通可达性布设采样点(图 1)，2018年8月对研究区进行采样，每个样地设立20 m×20 m样方，每个样方按照对角线采取3个重复样(包括环刀样2个)，采集样品深度为0 ~ 20 cm，共采集150个样品，每个样品由3个样点土样混合而成。土样去除植物根系和残渣后装入贴有标签的自封袋称重，带回实验室风干。使用马尔文(Malvern)激光粒度仪测量土壤中值粒径；使用碳氮分析仪(Elementar Vario Max)测定土壤全碳(TC)和全氮(TN)；环刀样主要用来测定土壤容重(BD)；有机质(SOM)测定采用重铬酸钾油浴加热法测定；电导率(EC)用EC快速测定仪测定；pH用pH计测定。

1.3 其他数据来源

气象数据：在中国气象科学数据共享服务网(http://cdc.cma.gov.cn)下载国家台站监测数据，通过Anusplin软件进行插值，得到研究区2018年年降水量(PER)与年平均气温(TEM)数据；归一化植被指数(NDVI)数据：本文使用的植被遥感数据MOD13Q1，来源于NASA网站(https://ladsweb.modaps.eosdis.nasa. gov/)，空间分辨率为250 m×250 m，时间分辨率为16 d，经过拼接、裁剪及最大合成法形成年NDVI数据；DEM数据：来源于地理空间数据云平台(http://www.gscloud.cn)，空间分辨率为90 m×90 m，经过计算得到海拔高度、坡度(Slope)、坡向(Aspect)、地形湿度(TWI)数据。

1.4 研究方法

本研究通过ArcGIS10.2软件提取点数据，利用Excel 2013软件对数据进行整理，利用平均值±3倍标准差去除异常点，通过GS+ Version 9软件进行地统计分析，对因变量进行插值，利用地理探测器模型进行影响因素分析。

1.4.1 地统计方法

1) 莫兰指数。莫兰指数反映研究区整体的空间模式，用于判断空间变量在区域内属于聚集分布还是孤立分布、研究属性值在空间位置上是否具有相关性^[17]。

$ I = \frac{n}{{\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {wij} } }} \cdot \frac{{\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {\left( {{x_i} - \overline x } \right)\left( {{x_j} - \overline x } \right)} } }}{{\sum {{{\left( {{x_i} - \overline x } \right)}^2}} }} $

(1)

式中：n为参与分析的土壤样品数目；x_i和x_j分别为样品i和j样品的观测值；$ \overline x $为所有样品观测值的平均值；w_ij为空间权重矩阵值。

2) 半方差函数。为探讨土壤全碳全氮的空间变异结构性和随机性，利用地统计学中的半方差函数分析其空间异质性变化。

$ \gamma \left( h \right) = \frac{1}{{2N\left( h \right)}}\sum\limits_{i = 1}^{N\left( h \right)} Z {\left[ {\left( {{x_i}} \right) - Z\left( {{x_i} + h} \right)} \right]^2} $

(2)

式中：γ(h)为半方差函数；N(h)为间隔距离等于h的样本点的对数；Z(x_i)为空间位置点x_i处指标的实测值；Z(x_i+h)为空间位置点x_i+h处指标的实测值^[18]。

3) 克里格插值。

$ Z({x_0}) = \sum\limits_{i = 1}^n {{\gamma _i}Z\left( {{x_i}} \right)} $

(3)

式中：Z(x₀)为研究区未知样本点的土壤全碳或全氮含量(g/kg)；Z(x_i)为研究区未知样点周围已知样点的土壤全碳或全氮含量(g/kg)；γ_i为第i个已知样点对未知样点的权重；n为已知样点的个数^[19]。

1.4.2 地理探测器

地理探测器是一种统计学方法，通过探测因变量空间分异性，揭示其背后驱动力及各因子的贡献率，能从庞大的空间数据库中提取有用的空间关联规则^[20-21]，由因子探测器、生态探测器、风险探测器及因子交互作用探测器4部分组成。

1) 因子探测器：用来分析影响因子的变化趋势对因变量的空间分异影响大小，用q值度量。其模型如下^[22]：

$ q = 1 - \frac{{\sum\limits_{h = 1}^L {{n_h}\sigma _h^2} }}{{n{\sigma ^2}}} = 1 - \frac{{{\text{SSW}}}}{{{\text{SST}}}} $

(4)

式中：q是度量空间分异性的指标；h=1, 2, …，L为分类数目；n_h和n分别为层h和全区的样本单元数，σ_h²和σ²分别为层h和全区的方差；SSW和SST分别为层内方差之和(within sum of square) 和全区总方差(total sum of squares)。q的取值范围为[0, 1]，q值越大，说明影响因子的解释力越强。

区域Y值方差的计算公式如下^[23]：

$ {\sigma ^2}{\text{ = }}\frac{1}{{{n_{h - 1}}}}{\sum\limits_{i = 1}^{{n_h}} {\left( {{Y_{h, i}} - \overline {{Y_h}} } \right)} ^2} $

(5)

2) 生态探测器：用于比较两因子对土壤全碳全氮空间分布的影响是否有显著差异，可判断两因子哪一个因子更具有影响力，以F统计来衡量^[24]：

$ q = \frac{{{N_{x1}} \times \left( {{N_{x2}} - 1} \right) \times {\text{SS}}{{\text{W}}_{x1}}}}{{{N_{x2}} \times \left( {{N_{x1}} - 1} \right) \times {\text{SS}}{{\text{W}}_{x2}}}} \leqslant $

(6)

$ {\text{SS}}{{\text{W}}_{x1}} = \sum\limits_{h = 1}^{L1} {{n_h}} \sigma _h^2, \;{\text{SS}}{{\text{W}}_{x2}} = \sum\limits_{h = 1}^{L2} {{n_h}} \sigma _h^2 $

(7)

式中：N_x1及N_x2分别表示2个因子的样本数量；SSW_x1和SSW_x2分别表示由2个因子形成分层的层内方差之和；L₁和L₂分别表示变量x1和x2数目。

3) 风险探测器：用于判断评价指标子区域间的属性均值是否有显著差别，用t统计量来检验^[25]。

$ t = \frac{{{{\mathop Y\limits^\_ }_{h = 1}} - {{\mathop Y\limits^ - }_{h = 2}}}}{{{{\left[ {\frac{{{\rm{var}} \left( {{Y_{h = 1}}} \right)}}{{{n_{h = 1}}}} + \frac{{{\rm{var}} \left( {{Y_{h = 2}}} \right)}}{{{n_{h = 2}}}}} \right]}^{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.} \!\lower0.7ex\hbox{$2$}}}}}} $

(8)

式中：Y_h为子区域h内植被NDVI属性均值；$ {n_h} $为子区域h内样本数量；$ \rm{var} $表示方差。

4) 因子交互作用探测器：用于识别指标之间的交互作用，即两因子之间共同作用(增加或减弱)或独立作用对土壤全碳全氮空间分布的解释力^[25]。

根据地理探测器数据要求，对自变量进行离散化，NDVI根据李晓香和陈媛丽^[26]的研究分为5个等级；Slope根据汤国安和宋佳^[27]的研究分为6个等级；海拔高度根据中国陆地高精度地貌类型的划分，分为3个等级；Aspect根据实际情况进行划分，分为10级；SOM划分为4个等级；其余因子均按照等间距进行划分，BD划分为5个等级、EC划分为7个等级、pH划分为3个等级、粒度划分为4个等级、TEM与TWI共划分为5个等级、PER划分为4个等级。各指标等级划分的具体结果见表 2。

2 结果分析 2.1 土壤全碳全氮含量描述性统计

从表 3中可以看出，TC含量的范围在8.38 ~ 112 g/kg，平均值为48.35 g/kg，根据第二次土壤普查养分分级标准，TC含量属于第一级；TN含量的范围在0.38 ~ 8.20 g/kg，平均值为4.00g/kg，根据第二次土壤普查养分分级标准，TN含量属于第一级。偏度及峰度是反映数据正态分布的两个指标，TN与TC含量的偏度均大于0，峰度均小于0，说明数据分布为右偏态，且较为平坦。从变异系数看，二者10% ≤CV≤100%，属于中等变异，K-S检验均在1% 的检验水平下，说明适合克里金插值。

2.2 土壤全碳全氮含量地统计分析

块金值与基台值之比表示随机部分引起的空间变异占系统总变异的比例^[28]。从表 4中可以看出，TN与TC含量的块金值/基台值比值在25% ~ 75%，属于中等变异，由结构性因素与随机性因素共同作用。

由图 2可知，TC与TN含量的空间自相关性随着自相关距离的增加逐渐减弱，TC含量的自相关距离在0 ~ 53 333.33 m，TN含量的自相关距离在0 ~ 70 366.67 m，其莫兰指数为正值，说明呈正相关，空间聚焦性较强；TC含量的自相关距离在53 333.33 ~ 106 666.67 m，TN含量的自相关距离在70 366.67 ~ 140 733.33 m时，其莫兰指数在0值上下浮动，说明空间自相关性较弱；TC含量的自相关距离在大于106 666.67 m，TN含量的自相关距离在大于140 733.33 m时，其莫兰指数小于0，说明TN与TC含量数据点存在空间孤立性。

从图 3中可以看出，TC与TN含量空间分布差异明显，高值区主要分布在研究区的南部，低值区主要分布在研究区的北部顶端，中值区主要分布在研究区的中部及南部，呈梯度变化分布，均一性较好，TN与TC含量呈现出东南向西北递减的趋势。

2.3 土壤全碳全氮含量的影响因素

利用因子探测对影响TC及TN含量的土壤属性、气象因子、植被因子、地形因子进行分析，可以看出，除了TWI及Slope因子的P值大于0.05，其余均小于0.05，说明除了这两个因子，TN及TC在其余因子分区的空间分布存在空间异质性且结果可信。自变量对TN含量的解释力大小排序为：SOM > BD > NDVI > EC > PER > Aspect > pH (DEM) > TEM > 粒度；自变量对TC含量的解释力大小排序为：SOM > BD > EC > NDVI > PER > Aspect > TEM (DEM) > pH (粒度)。可见，土壤内部属性对于TN与TC含量影响最大，其次为植被，再次为气象，最后为地形。

利用交互作用探测器识别影响因子之间的交互作用，从表 6中可以看出，任何两因素之间相互作用均大于单因素作用，且都为双因素增强；对于全氮，两因素交互作用解释程度大于0.6的因子对分别为SOM∩BD、SOM∩EC、SOM∩pH、SOM∩PER、SOM∩NDVI、SOM∩Aspect；对于全碳，有机质与其他所有因子的交互作用均在0.6以上。

利用生态探测器探测10个影响因素对土壤TN及TC含量解释程度的显著性差异，从表 7中可以看出，除了土壤内部属性对TN及TC含量的影响具有显著性差异，其他均无显著性差异。

风险探测用于判断影响因子子区域间的属性均值是否有显著的差别，用于探索土壤肥力最好的区域，据此，本研究得出，土壤TN及TC含量平均最大值出现区域为BD < 0.6 g/cm³，EC在0.3 ~ 0.35 µS/cm，pH在7 ~ 8，粒度在10 ~ 100 µm，SOM在100 ~ 150 g/kg，DEM在2 500 ~ 3 500 m(高中山)，PER在360 ~ 390 mm，TEM在–1 ~ 1 ℃，NDVI在0.95 ~ 1，TWI在6 ~ 9，Slope在15° ~ 25°，Aspect在北方。

3 讨论 3.1 土壤内部因素对土壤全碳全氮的影响

SOM与土壤全碳全氮关系较为紧密，相关性高达0.8左右；除SOM外，对土壤全碳全氮解释力最强的因子为BD，且与TN、TC呈负相关，BD通过影响土壤的粒度，进而影响土壤的含碳量，BD越大说明土壤孔度越小，土壤中含有的黏粒越多，土壤越紧实，不利于微生物活动，因此在BD < 0.6 g/cm³、土壤粒度10 ~ 100 µm时，TN、TC平均值最大；而EC可以间接地反映土壤含盐量，是植物和微生物活性的限制因素，影响土壤养分的转化，表征土壤酸化和次生盐渍化状况，同时也反映土壤的肥力状况。pH与TN、TC呈负相关，这可能是由于在酸性的环境中，微生物对有机碳的分解会受到抑制^[29]，研究区主要以中性或碱性土壤为主，因此，在pH 7 ~ 8范围内出现TN、TC平均值最大。

3.2 植被对土壤全碳全氮的影响

植被的组成、分布及群落结构都会造成土壤养分的变异^[30]，且地上生物量的大小直接影响土壤有机碳含量的高低，所以NDVI在0.95 ~ 1范围内出现TN与TC的平均最大值。研究区从南到北植被分布差异性较大，典型的植被类型有针叶林、阔叶林、针阔混交林及灌木林，而西北部地区海拔较高，主要以草甸、草原为主^[31]，中部主要以草原及灌丛为主，东南部主要以林地及草原为主，林地主要以祁连圆柏及云杉为主，分布于坡度较陡的山坡上，由于其植被分布的差异性，导致了土壤养分的带状分布。

3.3 气象因素对土壤全碳全氮的影响

根据土壤养分影响因素的各类研究报告^[32-33]可以看出，气候对土壤养分的影响占重要位置，尤其是表层土壤受气候因素的影响最大。本研究中，降雨对土壤养分的影响大于温度，研究区年降水量在340 ~ 400 mm，年蒸发量为1 496.3 mm，蒸发量是降水量的4倍左右。在强蒸发的作用下，水分成为了土壤养分变化的限制性因子，加大了土壤水分与土壤全碳全氮关系的紧密性^[34]，因此降雨与土壤养分的关系更为紧密。张华等^[24]得出在祁连山国家公园，年降雨量在320.7 ~ 385.8 mm、年均温在0.2 ~ 5.2℃时，最适宜植被生长，而本研究中年降雨量在360 ~ 390 mm、年均温在–1 ~ 1 ℃范围内TN与TC的平均值最大，两者研究结果相符，范围的差异是由于指标的等级划分不同造成的。

3.4 地形对土壤全碳全氮的影响

地形对于土壤的发育至关重要，能够影响水热条件的空间分配，改变区域内的气候条件^[35]。研究区主要以山地为主，山坡主要以林地为主，土壤较为肥沃，因此Slope在15°~ 25°范围内TN与TC平均值最大，表明在地形条件中坡向对于土壤养分的影响最大，且北方(阴坡)土壤TN与TC的平均值最大，这与袁杰^[15]对于祁连山南坡的研究结果一致，出现这一现象可能是由于阴坡蒸发小、气候阴冷潮湿，且土壤对水分的截留率较多，有较厚的枯枝落叶层，导致阴坡土壤相较阳坡土壤有机物质摄入来源多元化，相应的土壤肥力更好。

4 结论

1) 黑河源区土壤全碳与全氮含量的插值模型均为指数模型，插值的影响因素为结构性因素大于随机性因素，插值结果呈现出东南向西北递减的趋势，且呈梯度变化，这与黑河源区的植被地带性分布规律有密切关系，黑河源区林地主要分布在东南部而草甸主要分布在西北部与中部。

2) 通过地理探测器分析可以看出，SOM、BD、EC及NDVI是影响土壤全碳全氮含量变化的主控因素，可见，土壤内部因素的变化对于土壤全碳全氮变化起主要作用。

3) 土壤TN及TC平均最大值出现的区域土壤BD < 0.6 g/cm³，EC在0.3 ~ 0.35 µS/cm，pH在7 ~ 8，粒度在10 ~ 100 µm，SOM在100 ~ 150 g/kg，DEM在2 500 ~ 3 500 m(高中山)，PER在360 ~ 390 mm，TEM在–1 ~ 1 ℃，NDVI在0.95 ~ 1，TWI在6 ~ 9，Slope在15° ~ 25°，Aspect在北方。