土壤颗粒组成或机械组成是土壤重要的物理属性^[1-2]，是影响土壤肥力、作物产量^[3-5]、水分在土体内运移和储存的关键因素^[6]，同时也是陆面过程模型、水文模型和耦合陆面过程的大气模型的重要输入参数。

已有研究表明土壤颗粒组成受土壤有机质、气候、植被类型、地形地貌、土地利用等因素的影响^[7]。如罗绪强等^[8]研究发现喀斯特高原区贵州省清镇市土壤颗粒组成主要受土壤有机质的影响；朱丽东等^[9]研究发现洞庭湖地区土壤颗粒组成主要受气候因素的影响；夏江宝等^[10]研究发现红壤丘陵区土壤颗粒组成主要受植被类型的影响；庞龙辉等^[11]研究发现地形因素会影响土壤颗粒在空间上的分布；张宏和刘建军^[12]分析了在黄土沟壑区不同的土地利用方式下土壤颗粒组成的垂直分布特征。这些研究结果表明土壤颗粒组成受多种因素的影响，且存在着明显的区域差异。

近年来，机器学习方法被广泛应用于土壤颗粒组成的数字制图，如孙艳俊等^[13]采用基于人工神经网络的方法进行了土壤颗粒组成制图；刘超等^[14]采用支持向量机、模糊逻辑和决策树的方法对黑河流域土壤颗粒组成进行制图研究。随机森林作为一种新兴的机器学习方法，也开始被应用于数字土壤制图与环境因素之间的关系探究，如李璇等^[15]利用随机森林模型预测了祁连山区土壤砂粒含量与环境变量之间的关系；史静静等^[16]在黑龙江鹤山农场运用随机森林模型对影响包括土壤颗粒组成的土壤属性地形因子进行重要性排序并制图。

青藏高原是世界上海拔最高的高原，有“世界屋脊”之称^[17]，其气候具有独特的特征^[18]，是全球气候变化响应较敏感的地区^[19]，会对全球气候的变化产生重要的影响^[20]。此外，青藏高原分布着约为105 ~ 150万km²的冻土面积，其土壤具有一系列独特的特征^[21-22]，加之地形复杂，区域土壤颗粒剖面分布模式一般多是复杂多样，很难定量地预测^[15]，因此，目前对于青藏高原土壤颗粒组成的剖面分布模式及其与环境因素之间的关系的报道甚少。为此，本文选择青藏高原东部温泉地区，开展土壤调查并获取其土壤剖面数据，采用随机森林模型探讨各个土壤剖面颗粒组成的垂直分布模式及其与环境因素之间的关系，并希望得到影响其分布的主控因子，以期为青藏高原的地球系统模式科学研究和相关工程建设提供基础数据和科学依据。

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

青藏高原温泉区(35°06′ ~ 35°42′N，99°06′ ~ 99°42′E)位于青海省玛多县、兴海县、都兰县及玛沁县的交界处(图 1)，总面积约3 800 km²，其海拔梯度大，介于3 405 ~ 5 294 m，地形复杂多样。研究表明该区气温变化大，全年气温较低^[19]，年均气温-3.2 ℃，最高气温为18.5 ℃，最低气温为-32.5℃^[23]；年均降水量500 ~ 600 mm，属寒温带大陆性气候^[24]。根据野外实地调查，研究区以草原生态系统为主，研究区内以牧业为主，成土母质类型相对较少，以坡积物、冲积-洪积物为主，母岩有碳酸盐岩类、碎屑岩类、火成岩及变质岩类^[25]。

1.2 土壤剖面的布设与调查采样

由于研究区处于高山地带，受交通不便和地理条件的限制，数据采集不能按照简单的随机取样法进行均匀采样^[26]。温泉区样点的采集遵循先沿着G214国道进行，再根据沿线海拔、植被分布状况和地貌等因素进行局部调整的原则进行，这种基于专家知识和野外实际情况结合的采样方案在很大程度上可以减弱数据采集不均匀所造成的问题^[22]。根据采样方案共获得了73个样点及其详细的剖面、植被、地貌等信息。采用基于沉降法原理的SEDIMAT 4-12土壤粒径分析系统在实验室测定各个剖面不同层次土壤颗粒组成的含量及理化性质，粒度分级采用美国制，将土壤颗粒按大小分为黏粒(< 0.002 mm)、粉粒(0.05 ~ 0.002 mm)、砂粒(2 ~ 0.05 mm)3类^[27]。按照《中国土壤系统分类检索(第三版)》分类标准^[28]，该区分布有干润雏形土、潮湿雏形土、寒冻雏形土、永冻潜育土、干润均腐土5类土壤(亚纲)。

1.3 环境变量数据及其来源

土壤是气候、生物、地形、母质和时间综合作用的产物^[29]，时间因素不易表达，主要通过气候及植被因子间接体现^[30]。因此本文最终选择的环境变量主要包括气候、生物、地形和母质4个方面，其中气候因子主要选用了地表温度，考虑到青藏高原海拔高、气温变化大，根据汤懋苍等^[31]对青藏高原冷暖季的划分标准(9月至4月为冷季，5月至8月为暖季)，将地表温度分解为分季变量(冷季和暖季)^[22]。生物因子主要选用了NDVI；地形因子有高程、坡度、坡向、平面曲率、坡面曲率和地形湿度指数；母质因子主要是成土母质；此外，还加入了土壤类型和土地利用因子。这些环境因素可分为两类：可量化的和不可量化的。前者包括：高程、坡度、坡向、平面曲率、坡面曲率、NDVI、地形湿度指数和地表温度，其具体特征如表 1所示；后者包括：土壤类型、成土母质、土地利用，其类型和剖面个数分布如表 2所示。

采样点的经纬度和海拔数据由手持式GPS仪获取，研究区的高程、坡度、坡向、平面曲率、坡面曲率和地形湿度指数等地形因子数据是基于CGIAR- CSI提供的30 m SRTM-DEM数据运用ArcGIS10.2和3dMapper软件生成；NDVI数据由NASA提供的2000—2009年分辨率为250 m、16 d合成的MODIS-NDVI数据生成；地表温度数据由NASA提供的2003—2008年分辨率为1 km的MODIS逐日地表温度数据生成；成土母质、土地利用和土壤类型数据是在野外调查过程中获得。

1.4 随机森林模型

随机森林模型是Breiman于2001年在决策树算法的基础上发展而来，相较于其他分类方法，随机森林具有预测精度高、运行速度快、可避免过拟合、用袋外估计(OOB error)进行无偏估计等优点^[32]。随机森林模型可以通过R语言中的Random Forest软件包实现。

1.5 数据处理与分析

数据处理主要包括土壤剖面数据的层次标准化处理和各种环境因子的提取。由于野外调查获得的73个样点中，有15个样点的剖面深度≤1 m，为确保分析结果的可靠性，本文只对其余58个剖面深度≥1 m的标准化样点进行剖面分布模式分析，此外，由于测定的土壤质地是基于发生层，不同剖面划分的层次数目和各层次的深度不一，因此需要标准化处理。本文根据大多数据土壤剖面的发生层次情况，将标准化层次的深度分为5层，依次是0 ~ 15、15 ~ 30、30 ~ 60、60 ~ 100和100 ~ 300 cm。对剖面深度≥1 m的58个标准化样点，按照美国制的不同粒级，在Microsoft Excel 2016中进行标准化处理，其计算方法如公式1所示。

$ \sum\limits_{i = 1}^n {{\rm{P}}{{\rm{C}}_i} \times {{\rm{H}}_{\rm{i}}}} $

(1)

式中：PC_t-b为标准化层次t-b的某一粒级的含量，t、b分别代表该层次的顶部深度和底部深度；n代表剖面和发生层数目，PC_i代表发生层i的某一粒级的含量，H_i代表i层的深度在标准化层次t-b中所占的比例。

利用Microsoft Excel 2016和IBM Statistics SPSS 24.0软件对标准化后的土壤样点颗粒含量进行统计分析，再借助Origin 8平台进行土壤颗粒剖面分布模式的制图。环境因子的提取主要借助于ArcGIS10.2、3dMapper和ENVI软件，对可定量化的环境变量，采用随机森林模型分析其与土壤剖面分布模式之间的关系，对不可量化的环境因子进行统计分析。利用varImpPlot()函数对影响黏粒含量剖面分布模式的环境变量进行重要性排序^[33-34]，其中，Mean decrease accuracy(MDA)(把一个变量的取值变为随机数，随机森林预测准确性的降低程度)和Mean decrease gini(MDG)(计算每个变量对分类树每个节点上观测值的异质性的影响，从而比较变量的重要性)值越大，表明自变量越重要。

2 结果与分析 2.1 土壤颗粒组成垂直分布模式

表 3是58个剖面经过标准化处理后的砂粒、黏粒、粉粒含量的统计信息，可以看出，砂粒含量最多，平均含量在810 g/kg左右，变异系数为8.3% ~ 11.8%；黏粒居中，平均含量在100 g/kg左右，变异系数为36.8% ~ 66.2%；粉粒最少，平均含量在90 g/kg左右，变异系数为36.5% ~ 50.7%。从各个层级土壤砂粒、粉粒和黏粒的峰度和偏度值发现，基本服从正态分布。黏粒的变异系数明显高于砂粒和粉粒，其剖面分布的特征应更为明显。为了更好地理解青藏高原温泉区土壤颗粒组成垂直分布特征，选择变异系数大的黏粒进行分析。

对标准化后的58个样点在5个剖面深度处的黏粒数据进行统计分析，发现其可分为4种类型：递减型、先增后减型、先减后增型和不规则型(图 2)，其中递减型的分布模式样点最多，有19个，占32.8%；随后依次是先减后增型、不规则型和先增后减型，分别有17、16和6个样点，分别占29.3%、27.6% 和10.3%。递减型、先增后减型和先减后增分布模式的标准差线图和各自的剖面分布模式保持一致，不规则分布模式由于各个标准层数据规律比较差，其标准差线图呈先增后减型分布(图 2)。

2.2 土壤黏粒剖面分布模式与环境因素之间的关系

对影响黏粒含量剖面分布模式的环境变量进行重要性排序，结果见图 3和表 4。从图 3可以看出，MDA和MDG两个指标在同一剖面类型中对同一环境因素的排序不一定相同，为保障筛选的准确性，将同时位于MDA和MDG两个指标前3位的环境因素作为主控因子。结合图 3中的这两个特征值进行综合分析，结果显示，对于递减型分布模式，坡度、坡向和冷季地温是其主控因子；先增后减型的主控因子为暖季和冷季地温；先减后增型的主控因子为高程、地形湿度指数和NDVI；不规则型的主控因子为地形湿度指数、平面曲率和高程。

从表 4分析可知，递减型的成土母质主要为坡积物，而先增后减型、先减后增型和不规则型的成土母质主要是坡积物和冲积-洪积物，这说明递减型主要发育在山坡上，而其他3类主要发育在山坡下部或底部。4种模式在土地利用和土壤类型方面具有相似性，即先减后增型的土地利用和土壤类型分别是高寒草原和干润雏形土，而其他3类均分别为高寒草甸草原和干润雏形土与干润均腐土。高寒草甸草原相对于高寒草原来说，其分布的位置往往水分条件更好；雏形土是剖面发育程度较低的未成熟土壤，在本区域内大面积分布的主要是因为温泉地区年均气温较低，且降雨相对较少，导致土壤发育程度较弱。而均腐土是指具有暗沃表层和均腐殖质特性的土壤，在本区域内往往发育在水分条件相对较好、高寒草甸分布广泛的地方。从两种土地利用类型和土壤类型对比可以看出，先减后增型剖面在本区域内主要分布在相对较干的位置。

2.3 黏粒含量与不同剖面分布模式的典型案例

结合野外考察和实验室对各个层次土壤理化性质的分析发现，递减型剖面分布模式的主要特征是土体上部以较细的黏粒和粉粒为主，下部以较粗的寒冻风化物(块状砾石)为主，黏土的含量随着深度的增加呈现递减的趋势。如图 4A为该模式的典型剖面，其位于山体阴坡中部，其细颗粒物质在重力和较高雨水共同作用下，易向外置更低的地区迁移，在垂直分布上黏粒含量就表现为随着深度的增加而减少。这也进一步映证了坡度和坡向作为递减型模式的主控因子，通过影响土体的稳定性和地表水流的聚排^[35]对土壤颗粒的分布产生了影响；而冷季地温过低会抑制植被的生长，进而影响研究区土壤的发育，递减型土壤剖面上部土层较薄，下部基本上是由寒冻风化物组成的坡积物。

先增后减型和先减后增型剖面分布模式主要特征是由于受到坡积、洪积作用的影响，土体来源较复杂，导致其剖面分布模式也复杂多样。如图 4B为先增后减型的典型剖面，位于山坡底部，土体由来自不同时期的堆积物构成，土体较厚，土体中上部和下部黏粒含量较少，中间黏粒较多，下面发育多年冻土，总体上土壤发育较好，地表温度(暖季和冷季)是其主控因子。图 4C为先减后增型的典型剖面，其所在的位置海拔更低，高程与地形湿度指数反映出其是个低洼的汇水区域，NDVI反映了其植被状况较差；土体上部以块状的细壤土为主，中部以砂壤土为主，下部以单颗粒的砂和黏土为主；剖面中植被根系较少，底部有积水，因下面没有发育多年冻土(隔水层)，土壤持水能力差而导致植被较少，这与前面得到高程、地形湿度指数和NDVI为其主控因子一致。

不规则型分布模式的母质也为坡积和冲积-洪积物，剖面主要表现为土-石交错排列的特征。如图 4D是不规则型的典型剖面，其位于海拔4 250 m的坡积裙底部，经过较长时间的坡积和洪积相互作用，大颗粒的砂石和细颗粒的土体就交错分布；此外，图 4D剖面下部有冻融扰动现象，判定其下面发育多年冻土，根据前人^[36]的研究，温泉地区的多年冻土下界高程约为4 090 m。前文得到不规则型分布模式的主控因子主要为地形湿度指数、平面曲率和高程，地形湿度指数能有效地指示土壤内水分运动方向，平面曲率影响着物质的加速集中或减速分散，而高程则可指示多年冻土发育状况。

3 讨论

一般而言，土壤质地剖面垂直分布模式受成土母质、气候、植被类型、地形地貌、土地利用等多种环境因素的影响。对于地势平坦、母质连续、区域气候变化不大的地区，其土壤质地的剖面分布模式较为简单，多趋向于均一型；李璇等^[15]认为对于地形复杂、母质不连续的区域，其颗粒组成的剖面分布模式也往往比较复杂。但本文研究结果进一步显示，在母质相对连续单一的高山高寒地区，受海拔梯度变化大、年平均温度低和冻融循环剧烈等作用影响，其土壤质地剖面分布模式也会复杂多样。

温泉地区的成土母质主要为坡积物和冲积-洪积物，但分布着4种类型的土壤质地剖面模式，结合前面得到的主控因子来看，除了受到成土母质的影响外，其中地形和气候因素(地温)也对剖面模式的分布起着决定性的作用。先看地形因子，其中坡度和坡向是递减型模式的主控因子，主要是因为温泉地区有一系列相间排列的山地和断陷盆地，致使区内地形起伏大，山坡上部的细颗粒物质会在重力作用下向下方移动，覆盖在山坡中部的寒冻风化物上，从而形成了坡积作用，这一点也可以从递减型的成土母质为坡积物得到验证；此外，高程和地形湿度指数均为先减后增型和不规则型的主控因子，可以看出这两种模式分布位置都在山坡底部，都是汇水区域，区别在于不规则型分布的海拔位置更低，位于季节冻土区，而先减后增型位于多年冻土区，这主要是因为温泉地区海拔梯度大(3 405 ~ 5 294 m)，地形复杂多样；再看气候因素，平均地温均是递减型和先增后减型的主控因子，足可以看出气候因素对本区域土壤质地剖面模式的影响，主要是因为温泉区处于高寒地区，年平均气温为-3.2 ℃^[20]，且发育着大面积的多年冻土，据前人研究^[36]，温泉地区多年冻土的面积约占整个研究区总面积的76%，低温和剧烈的冻融循环都对本区域土壤的发育产生了较大的影响，也对本区域的土壤质地剖面模式的分布产生了影响。

相比于前人在其他区域研究得到的相似结论，如朱丽东等^[9]研究发现洞庭湖地区土壤颗粒组成主要受气候因素的影响；卢宏亮等^[37]研究发现安徽省土壤黏粒含量在空间的分布主要受地形因素的影响；庞龙辉等^[11]研究发现地表温度和地形因素会影响青海省土壤颗粒在空间上的分布。可以发现影响洞庭湖地区土壤质地的气候因素主要是指该区域的季风风向和大小，因为该区域的土壤颗粒组成具有较强的风成特性；影响安徽省土壤质地的地形因素主要是高程、山谷平坦指数和平面曲率，是因为该区域内地形地貌复杂，既有平原和丘陵，又有大别山这样的高山区；而影响青海省土壤质地的环境因子主要是地表温度和高程，本文研究结论与此一致，因为本文研究区也位于青海省，属于青藏高原范围内，海拔梯度大、地形复杂多样、年均温度低，导致了气候和地形因素是其主控因子。

综合以上分析，可以看出土壤颗粒组成受多种环境因素的影响，且存在着明显的区域差异。而青藏高原作为我国的高山和高寒区，其主控因子既有和其他地区相似的地形因子，也有自己独特的气候因子。整体来说，地形和气候因子是共同影响青藏高原温泉地区土壤黏粒剖面分布模式的决定性因素，这也是青藏高原所独有的特征。

4 结论

青藏高原温泉地区土壤砂粒含量最大，但黏粒的变异系数最大。黏粒的剖面分布模式可分为递减型、先增后减型、先减后增型和不规则型4种类型。

黏粒递减型分布模式的主控因子是坡度、坡向和冷季地温，先增后减型分布模式的主控因子是暖季和冷季地温，先减后增型分布模式的主控因子是高程、地形湿度指数和NDVI，不规则型分布模式的主控因子是地形湿度指数、平面曲率和高程。气候和地形因素是影响土壤黏粒剖面分布模式的决定性因素。