土壤养分是土壤提供的植物生长所必需的营养元素^[1]，是土壤生产力的基础，了解土壤养分的空间变异，对于研究土壤退化、生物多样性缺失、气候变化和制定土地管理计划具有重要意义^[2-3]。我国西部省区，尤其是高寒山区，土壤景观关系复杂，土壤空间变异程度高^[4]，但缺乏准确的土壤空间信息，制约着农牧业的科学生产与管理。

数字土壤制图是将土壤-景观模型通过计算机辅助生成的以栅格数据为表达形式的制图方法，可反映较为准确的土壤空间信息。制图模型主要包括地统计、专家知识和机器学习方法，研究内容集中于土壤的多种属性。在区域从小到大、景观从简单到复杂、制图结果从二维到三维(多维)的制图研究中，环境变量选取的是相对较少且易获取的协变量集，主要集中于地形因子与气候因子以及用遥感影像来表征的生物因子。涉及大面积区域时，成土母质这一因子的获取较难，往往被忽略；同时，表征土壤和环境之间关系的一些成土因子定量化很难表示(如时间因子)^[5-12]。因此，要实现大尺度、高精度的制图，开发挖掘丰富的环境协变量数据库是一个重要的突破口。在现有的机器学习(集成)模型甚至是人工智能模型中，只能反映小区域、局部景观环境与土壤的关系，模型对于与研究区环境条件相似的地区可能适用性较好，对于复杂的大面积区域适用性较低^[12-14]。要精确地模拟土壤性质和过程，有必要开发具有处理不同区域土壤与环境之间非线性关系的针对性的模型。目前国内数字土壤制图工作主要集中在我国中部和东部地区，地形梯度相对较小，可达性较好，土壤调查样本通常数量较多且分布相对均衡。例如，杨煜岑等^[5]在陕西省蓝田县用多元线性回归模型预测了陕西省蓝田县土壤养分的空间分布；Dong等^[6]将机器学习算法与基于地块的DSM框架相结合，对宁夏冲积平原农业区表层土壤养分(土壤有机质、全氮等)特征进行了预测；Liu等^[7]用地表动态反馈方法预测了江苏省中部平原地区土壤有机质含量等属性的空间分布。然而，在我国西部大面积高寒山区，路网不发达，可达性较差，土壤调查样本较为稀疏，已有工作主要集中在景观或流域尺度以及单一土壤养分属性(如有机碳或全氮)，如代子俊等^[8]针对湟水流域分析了近30 a的土壤全氮的时空变异影响因素；张欢和高小红^[9]基于湟水流域土壤有机质属性对比分析了地统计学进行预测制图的精度；周伟等^[10]在三江源地区将原始光谱反射率数据及其不同数据变换形式下的光谱分别与土壤有机质含量进行了相关分析；庞龙辉等^[11]对青海省土壤全氮、有机碳、粉粒和pH进行了预测制图。如何对大面积高寒山区土壤养分空间分布准确地预测，并在此基础上面向农牧业进行合理的生产管理分区，是一个值得研究的问题。

本文以地处青藏高原东北部的青海全省为研究区，将机器学习方法与成土环境协同变量相结合，探讨土壤养分(全氮(TN)、全钾(TK)、全磷(TP))含量的空间预测制图及生产管理分区，以期为青海省土壤养分空间管理及农牧业可持续发展提供必要参考。

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

青海省地处我国西北部内陆腹地，地理位置为31°39′ ~ 39°19′N、89°35′ ~ 103°04 E′，总面积72.23万km²，大部分区域为高寒山区，地势总趋势为南高北低^[15]，海拔介于1 686 ~ 6 477 m，大部分地区海拔在3 000 m以上。地貌上复杂多样兼有高原、山地、草原、湖泊；气候上，虽然总体上属于高原大陆性气候，但属暖温带向亚热带过渡型，不同海拔高度和不同地区气候有明显差异，年均温度–5.1~ 9.0 ℃，年均降水15 ~ 750 mm，降水量由东向西逐渐递减。按大地构造全省大体可分为三类不同地形区：南部青南高原、西部柴达木盆地和北部与东部平行岭谷区^[4]。区域内成土母质类型较多，成因各异。受地形、气候、土壤环境的影响，草地面积占全省的3/5，畜牧业得以发展，全省的森林、草原、湿地等在应对气候变化上有显著的作用，是我国重要的生态屏障。研究区采样点空间分布见图 1。

1.2 数据来源与处理

土壤剖面数据来源于科技部基础专项“我国土系调查与《中国土系志—青海卷》编制”项目中的205个样点，主要包括采样点位置和景观信息以及土壤理化性质(TN，TP，TK，SOC，pH，黏粒、粉粒、砂粒含量等)。土壤TN采用重铬酸钾、硫酸消化–蒸馏法(凯氏蒸馏法)测定；TK使用碱熔–火焰光度法测定；TP使用碱熔–钼锑抗比色法测定^[16]。对剖面各发生层测定的数据进行深度加权，生成表层(0 ~ 20 cm)的TN、TK、TP数据。

环境变量的选取是以土壤发生学理论为依据，主要依照土壤形成的五大要素(母质、地形、生物、气候和时间)进行选取代表性环境变量。由于目前尚无统一的定量指标来描述时间信息，因此本研究暂不考虑时间类因子。用土壤类型(soil types)近似表征母质的差异；地形因子在土壤发生与演变过程中影响着地表物质能量的循环转换；遥感数据具有实时、易获取的特点，其特征波段及不同的组合可以用来表征生物景观信息的综合环境(地表景观信息)。气候因子(降水量、气温)是直接影响成土过程中的水热条件。本研究通过表层土壤属性及其与环境变量之间的相关性确定了比较关键的21个环境变量。其中，地形因子：DEM数据来自于资源环境科学与数据中心(http://www.resdc.cn/)，基于DEM数据利用SAGA软件(https://sourceforge.net/)计算提取海拔(Elev)、坡度(Slp)、坡向(Asp)、剖面曲率(Curprf)、地形湿度指数(Twi)、地形开放度(Openness)、多尺度山谷指数(Mrvbf)；气候数据：年均温(MAT)、年均降水量(MAP)、气温季节性(Ts)、昼夜温差(Mdr)、潜在太阳辐射(Radition)数据来自于世界气象数据库网^[17] (http://worldclim.org/version2)；土壤类型(Soiltype)和土地覆盖类型(Lucc)均来源于资源环境科学与数据中心(http://www.resdc.cn/)；植被指数均值(Evm)，植被指数标准差(Evs)，白天地表温度(LSTday)，夜间地表温度(LSTnight)，7、8月地表温度(LST78)和地表反射率(band4/7)的处理是在Google Earth Engine (GEE)云平台通过JavaScript编程完成，选取时间为2015年7—8月的MODIS影像数据集，处理包括影像去云、掩膜裁剪、投影变换。所有栅格环境变量均重采样(最近邻法)为1 km。

1.3 土壤养分空间预测方法

随机森林^[18]是通过建立多个决策树并将其融合成一个更加准确和稳定的模型，能够根据样本之间的相似性来衡量特征交互，广泛应用于分类和回归。随机森林建模在R.4.1中Random Forest包进行^[19]，该模型主要有3个参数需要定义，即树的数量(ntree)、叶节点包含最少样本的个数(nodesize)、随机选取变量的个数(mtry)。使用caret包对nodesize和mtry两个参数进行优化，优化值分别是5和3，ntree采用默认值500，足以产生稳定的预测。%IncMSE作为判断预测变量重要性的指标，该值越大表示该变量的重要性越大，即对土壤养分含量的影响度就越大。

本文中制图预测精度采用留一交叉验证法LOOCV(leave one out cross validation)。其原理是：在个数为N数据集中随机抽取1个样本去验证N-1个数据所训练的模型，重复N次取平均值^[20]。评价指标包括：决定系数R²、一致性指数(CC)、均方根误差(RMSE)、预测偏差值(ME)。其中，R²表示模型对自变量解释的比例，取值范围0 ∼ 1；CC表示预测误差的方差值占实测值的比，取值范围0 ∼ 1，1表示预测结果最佳，0表示结果不太符合；RMSE是预测值与实测值之间的均方根误差；ME是预测值与实测值之间误差值的平均，反映预测结果情况是否存在高估或低估。公式如下：

$ {R^{\text{2}}}\frac{{\sum\limits_{i = 1}^N {{{\left\{ {p(xi) - \bar p} \right\}}^2}} }}{{\sum\limits_{i = 1}^N {{{\left\{ {v({x_i}) - \bar p} \right\}}^2}} }} $

(1)

$ {\text{CC}} = \frac{{2 \cdot R \cdot {\sigma _p} \cdot {\sigma _v}}}{{{\sigma _p}^2 + {\sigma _v}^2 + {{\left( {\bar v - \bar P} \right)}^2}}} $

(2)

$ {\text{ME}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left\{ {p({x_i}) - v({x_i})} \right\}} $

(3)

$ {\text{RMSE = }}\sqrt {\frac{{\text{1}}}{N}\sum\limits_{i = 1}^N {{{\left\{ {p({x_i}) - v({x_i})} \right\}}^2}} } $

(4)

式中：N为样本的数量，$ p({x_i}) $为样本的实测值，$ v({x_i}) $为样本的预测值，$ \bar P $为样本实测值的平均值，$ \bar v $为预测值的平均值，$ {\sigma _p}^2 $为样本实测值的方差，$ {\sigma _v}^2 $为样本预测值的方差。

1.4 土壤养分分区评价方法

投影寻踪是将高维数据向低维空间投影，通过分析低维空间的投影特性来研究高维数据的特征，找到反映数据结构特征的最优投影，是处理多指标问题的统计方法^[21]。其具体思路是将影响问题的多因素指标通过投影寻踪分析得到反映其综合指标特性的最优投影特征值，然后建立投影特征值与因变量的一一对应关系，通过分析一维的投影值来对样本做出更加合理的分级和评价^[21]。

首先，假定$\left\{ {{{{x_{ij}}} \mathord{\left/ {\vphantom {{{x_{ij}}} {i = 1,\; \cdot \cdot \cdot ,\;n;\;j = 1,\; \cdot \cdot \cdot ,\;m}}} \right. } {i = 1,\; \cdot \cdot \cdot ,\;n;\;j = 1,\; \cdot \cdot \cdot ,\;m}}} \right\}$是样本集，${x_{ij}}$是第i个样本的第j个指标值，使用每个指标数据的最大值对该指标进行标准化处理，以计算$ \frac{{{x_{ij}}}}{{{x_{j\max }}}} $得到标准化后的样本数据$x_{ij}^*$。

然后，使用数量为x的多维单位向量a(a₁, a₂, …, a_m)表示某一投影方向，用公式(5)可计算样本i在该方向上的投影值，用公式(6)计算投影值的标准差S(a)，用公式(7)计算投影值的局部密度D(a)，用投影值的标准差S(a)和局部密度D(a)的乘积表示为投影指标函数Q(a)，即公式(8)，目的是当把m维数据综合成某一方向上的一维投影值时，投影值散布特征为局部投影点尽可能密集，而整体上各个点团之间又尽可能分散^[21]。

$ {z_i} = \sum\limits_{j = 1}^m {{a_j}x_{ij}^*,\;i = 1,2,3,\; \cdot \cdot \cdot ,\;n} $

(5)

$ S(a) = \sqrt {\frac{{\sum\limits_{i = 1}^n {{{({Z_i} - {{\bar Z}_i})}^2}} }}{{n - 1}}} $

(6)

$ D(a) = \sum\limits_{i = 1}^n {\sum\limits_{k = 1}^n {(R - {r_{ik}})f(R - {r_{ik}})} } $

(7)

式中：$\bar Z$为投影方向上投影值的平均值，${r_{ik}} = \left| {{z_i} - {z_k}} \right|$，R为局部密度的窗口半径，可取$ {r_{\max }} + \frac{m}{2} \leqslant R \leqslant 2m $，函数f为单位跃阶函数，R–r_ik≥0时，其值为1，R–r_ik < 0时，其值为0。

$ Q(a)=S(a)D(a) $

(8)

考虑到Q(a)值与投影方向a有关，当Q(a)取最大值时的a方向是最能反映数据结构特征的方向，因而利用遗传算法，通过求解投影指标函数优化问题识别最佳投影方向，公式如下

$ \left\{ {\begin{array}{*{20}{c}} {\max Q(a)} \\ {s.t.//a// = 1} \end{array}} \right. $

(9)

最后，基于最佳投影方向采用公式(5)计算全国土壤养分含量分级标准值的投影值，建立基于投影值的等级评价和分区方法。依据全国第二次土壤普查的分级标准(表 1)^[22]，将土壤TN、TP、TK的预测结果进行叠加分析、归一化处理，基于DPS9.5软件(http://www.dpsw.cn/)计算最优投影方向以及最优投影方向的投影值，较为客观地对土壤养分分区进行划分。

2 结果与分析 2.1 表层土壤养分的统计特征

表 2为表层样本的TN、TP、TK的描述性统计情况。Kolmogorov-Smirnov检验结果显示，TN、TP经对数转换后均符合正态分布，TK含量符合正态分布。样点表层土壤TN含量范围0.04 ~ 13.50 g/kg，平均值2.56 g/kg；TP含量范围0.49 ~ 3.39 g/kg，平均值1.52 g/kg；TK范围4.67 ~ 26.95 g/kg，平均值17.23 g/kg。从变异系数来看，表层土壤养分的变异系数均处于0.1 ~ 1，表明土壤养分均呈中等程度变异，其排序为TN > TP > TK，其中TN的变异程度最大，达到0.94。

2.2 土壤养分预测制图与结果

利用随机森林模型预测青海省表层土壤养分空间分布(图 2) 表明，土壤TN预测值的范围为0.20 ~ 8.67 g/kg，均值2.48 g/kg，在空间上整体呈现东南向西北递减的趋势，其中青海省东南、中部和东北部土壤TN含量较高，西北部及湟水谷地的含量较低。土壤TP的预测结果范围0.75 ~ 2.41 g/kg，均值1.33 g/kg，空间分布上的趋势是东高西低，土壤TP含量高的区域是在玉树地区的东南部、东北部的祁连山地区和全省东部，而西南部可可西里地区、西北部的茫崖地区和中部昆仑山区含量较低。土壤TK预测的结果范围12.70 ~ 22.25 g/kg，均值16.94 g/kg，整体呈现分布趋势为东南高西北低，空间分布上较为分散，土壤TK含量高值主要分布在青南高原和祁连山山区，在唐古拉山中段地区也呈现高值，低值主要分布柴达木地区、昆仑山脉南部、黄南州的西北部和青海省西南部的格尔木市地区。

交叉验证结果(表 3)表明，土壤TN、TP、TK模型拟合度R²比较高，分别为0.89、0.85、0.82，即可以解释土壤养分含量的80% 以上；从一致性系数(CC)来看，三者均在0.87以上；预测偏差(ME)三者接近0值，说明预测的值有被低估。由此，机器学习方法与成土环境协同变量相结合对青海省表层土壤养分的预测，可以很好地解释土壤养分的空间变异。

2.3 土壤养分空间变异的影响因素

图 3显示了土壤养分空间预测中环境协同变量的相对重要性排序。对于青海省表层土壤TN的空间分布格局，植被指数均值(Evm)为主导因子，年均降水量(MAP)次之。通常在空间范围较大研究区域内，气候因素的地带性可以直观反映作物的生长状况，间接影响着土壤养分元素的迁移转化过程。在青海省东北部的祁连山地、湟水谷地及果洛玉树高原区域的降水量较高，土壤水分条件好，植被覆盖度高，且随海拔增加，覆盖率增加，地表的枯枝落叶增加，腐殖质层较厚，土壤微生物活性较强，土壤养分来源更加丰富，积累的TN含量普遍较高。在青海省西部地区，气候干旱炎热，降水量稀少且地表温度高，水热条件的不均衡导致植被覆盖度小，影响着土壤中微生物的活动，TN含量低。其次，西北区域分布着广泛大面积的盐渍土、风沙土、砂质土等，且其土壤保水保肥能力较差，土壤养分分解速度快，植物受高温干旱和盐分制约生存，氮素的来源和累积弱，土壤TN较低。

从TP来看，植被指数均值(Evm)、地表覆盖(Lucc)和海拔(Elev)对其含量重要性排名前三。气候、地形地貌等环境因子对土壤中磷元素的循环产生重要影响，玉树地区的东南部、东北部的祁连山地区为青海省的植被覆盖度高的区域，土壤中有机质高，可减少磷的吸附位点，增加磷的移动性，其次土壤磷形态和含量受土壤微生物活性的直接影响，土壤中的微生物，能够通过其生理活动促进土壤中含磷化合物中的磷释放，从而改变土壤中磷元素形态和含量。中北部柴达木盆地，存在绿洲所以出现高值；西南部可可西里地区、西北部的茫崖地区和中部昆仑山区，气候干旱，同时受海拔、风力的作用增强，土壤中砂粒含量高，温度流失较快，土壤磷素的转换率低。青南高原河谷地区降水量相对充足，此外海拔(Elev)对降水(MAP)分布模式产生影响，进而土壤pH空间分布也受之影响，同时气温受海拔的影响，在此区域的蒸发量远低于降水量，使得土壤淋溶过程加速，促进土壤酸性升高，在中、酸性的土壤中，磷素形态转换速率快，土壤TP含量较高；西部柴达木盆地气候干燥，蒸发量大，容易形成盐沼，导致碳酸盐集聚，土壤TP含量偏低；在东北部的河湟谷地区域土壤TP含量高，局部变异也大，该区域是青海省主要农业区，长时间的人类农耕活动作用，导致土壤碱性增强，土壤TP含量高。

对于TK而言，植被指数均值(Evm)、年均降水量(MAP)、昼夜温差(Mdr)的影响作用很大。青海省南部的玉树、东南部的果洛地区和东北部的祁连山山区植被覆盖度高，植被覆盖较好，土壤受淋溶作用相对较弱，加之林地土壤中的充足草木枯萎被微生物分解，有利于土壤速效钾的涵养与保持。中部的昆仑山区和唐古拉山中段地区的昼夜温差(Mdr)较大，风雪气候常见，存在严重的土壤侵蚀，由于土壤全钾主要是由矿物钾组成，土壤钾素会随土壤侵蚀流失，土壤中TK含量较低；在西南部的格尔木地区，海拔较高，受风力与水流的风化作用较强，风化物中又富含较多的云母、长石等含钾结晶体，富钾风化残积物混入土壤后使得含量变高。潜在太阳辐射(Radiation)、年均降水量(MAP)等协变量因子对土壤的湿润和淹水时间持续的长短起到关键的作用，进而影响钾素固定和释放过程，土壤失水变干，会产生释钾现象，青海省地形起伏度大，同时受气候和人类活动影响，由此全省的TK含量局部上有所变异。

2.4 土壤养分分区结果

从表 4可知，土壤养分分级标准值对应的投影值分别为0.141、0.234、0.411、0.682和0.814，据此对土壤养分投影值划分等级，其结果为：1级(< 0.141)、2级(0.141 ~ 0.234)、3级(0.234 ~ 0.411)、4级(0.411 ~ 0.682)、5级(0.682 ~ 0.814)、6级(> 0.814)。

图 4为青海省土壤养分分区图，表 5是每个等级和分区所占面积的比重。土壤养分等级划分为极低、低、中、中上、高、极高6个等级，管理分区分为6个分区(Q1 ~ Q6)。全省土壤养分等级划分的大致空间格局为：整体上呈现东南-西北均匀递减趋势，极高的主要分布在玉树的南部地区、果洛地区和东北部的祁连山，其分布区域植被覆盖度高；昆仑山段也存在极高的分布，面积比例占全省的4.18%，这表明该地区的土壤人为活动较少，主要是自然发育的土壤，主要受气候因子和地形因子影响作用较大；养分等级高的区域面积比例占全省的面积的13%，其空间分布不均，主要分布在海北州、格木市飞地、玉树州的西南部和黄南州的西北部等地区，这些地区有大面积的林区、草场和耕地，同时受人类活动作用影响比较大；养分等级中至中上的区域主要分布西北柴达木盆地周边、东部湟水谷地、青海湖周围以及玉树州的西南区域，占面积比例65%，这与该区域内的地形坡度变化及人工施肥量有关；土壤养分低等级的区域主要分布在柴达木盆地、西南的可可西里区域、扎陵湖和鄂陵湖周边和海南州的中北部，占面积比16%。土壤养分等级从高到低的面积比例为：4︰13︰40︰24︰16︰3，养分等级低和极低的总面积占全省的19%，表明青海省土壤养分属于中等级及偏上，主要分布在全省的中南部和东北部。因此，需要因地制宜，科学地划分土壤利用区域，调整放牧和草场施肥结构，加强对土壤养分含量高的区域的管理与调整，恢复土壤养分含量低的区域，增加生态保护措施。

3 讨论

通常在研究区面积相对较小、景观相对单一、样点较多且分布相对均衡的条件下，只利用样点而不用环境变量进行简单插值预测，预测效果比较好。如张铁婵等^[23]和石小华等^[24]采用空间插值法分别对榆林市榆阳区和周至县北部猕猴桃适生区的土壤养分元素进行空间预测，在样点布设大的地区，插值精度很好，反之插值精度不是很理想。插值结果受样点布设密度的影响较大，而预测精度受不同的插值方法的影响较小。对于大面积高寒地区，样点稀少且分布不均，显然单纯基于样点的插值不适合，需要发挥环境因素辅助变量的作用。本文研究发现，通过把丰富的成土环境因素变量与复杂的机器学习方法结合，能够在大面积高寒山区实现较为准确的土壤养分预测制图，捕捉土壤空间变异。植被能有效地指示土壤养分在大面积区域的空间变异。植被指数均值(Evm)对青海省表层土壤养分的空间变异具有主导作用。Pouladi等^[25]在丹麦日德兰半岛，利用NDVI、SAVI、DVI、RVI和地形相关变量预测表层SOC含量，发现这些遥感变量解释了89% 的SOC空间变化。Wang等^[26]在辽宁省表层土壤含量制图中发现结合仅利用遥感数据的随机森林模型预测效果明显好于仅使用地形和气候变量的随机森林模型。在后续工作中，在景观复杂且草地大面积发育的青海省，遥感数据作为主要环境变量的土壤制图值得探究。地表温度是地表对于太阳辐射的响应或反馈^[27]，可以作为环境变量有效地对土壤环境的空间差异进行揭示，对于本研究的高寒山区，地形地貌复杂多变，降水少，植被以牧草为主，这种反馈将主要取决于太阳辐射、地形、降水、植被和土壤的综合作用。这也说明该区域地表温度(LST)和年均降水量(MAP)是影响表层土壤养分空间变异的主要环境因子。地形因子的贡献率在青海省表层土壤养分的空间变异中相对较低，主要的原因可能是在景观和小流域尺度上，母质、气候和植被相对一致，地形通过径流、土壤侵蚀间接影响土壤属性的空间分布^[28]，地形是影响小面积区域土壤属性的主要因素，如孙孝林等^[29]和杨琳等^[30]在小区域土壤养分预测中只考虑地形属性变量就可以进行有效的土壤制图。在大面积高寒山区，土壤属性受到多种环境因素的协同影响，如母质、气候、地形和植被等。各成土因素变量对土壤养分都有贡献，有较强的成土因素综合作用。土壤属性与环境变量之间存在非常复杂的关系，这与小面积区域土壤制图不同。在后续的工作中，需要发掘更为有效的成土环境因素协同变量，来提高高寒山区数字土壤制图的准确性。

另一方面，预测模型选取应考虑具体的土壤样点数据和土壤景观情况。大面积高寒山区地貌景观复杂多样，土壤空间变异程度高，样点比较稀少且空间不均衡，导致这类地区数字土壤制图结果难度往往比较大，简单景观是指某个成土因素主导土壤空间分布，土壤-环境关系近于线性。复杂景观是指多个成土环境因素，土壤-环境关系是非线性的、非平稳的^[31]。对于前者，线性回归可以较好地解释协变量(如地形属性)与土壤属性之间的关系。对于后者，可以使用机器学习算法来模拟土壤变化。青海省作为大面积高寒地区，土壤-环境关系复杂，因此预测模型的选取至关重要，张欢和高小红^[9]比较了地统计中协同克里格、地理加权回归克里格以及结合土地利用类型的均值、中值修正协同克里格等方法对湟水流域土壤有机质属性的预测精度，忽略了复杂的土壤环境关系，在空间分布上能体现其变化规律，但预测值结果存在预高估和低估现象甚至负值。地统计学模型在土壤样本均匀分布且满足空间自相关的条件下，模型的预测精度较好，这在地形复杂、土壤调查样本较为稀疏的青海省并不适用。

土壤的利用和管理需要对土壤养分进行生产管理分区。通常土壤属性指标和分区方法的选取不同，会忽略土壤综合养分的空间差异、土壤质量因子的变异性以及土壤环境的特殊性，进而管理分区结果也因此产生不同，大部分学者采用了主成分分析、地统计、神经网络和模糊聚类模型等方法结合土壤属性的单因子来进行分区管理，黄婷等^[32]通过主成分分析来计算土壤综合质量指数进行了线性回归分析，忽略了土壤-景观关系的复杂多变性；朱昌达等^[33]运用地统计方法分析土壤属性空间变异的同时采用模糊聚类来进行分区，未能考虑土壤属性的综合作用，计算分区指标权重和综合值。由于土壤综合养分具有多维性和综合性，而综合定量化研究在较大程度上可以避免研究者主观因素的影响，已然成为土壤综合养分分区与管理的一个趋势^[34-36]。对于大面积高寒山区，土壤养分与环境变量之间的关系较为复杂，需要建立土壤养分与评价值之间较为准确的关系，重视土壤养分的综合指标的选取信息，来提高管理分区结果的准确性，本研究基于数字土壤制图得到的土壤养分分布图，采用投影寻踪法对土壤养分进行了管理分区。尽管投影寻踪法减少了人为因素干扰，提高了养分管理分区的客观和准确性，但是不同投影方向优选算法可能会影响评价结果，有待进一步探讨。

4 结论

1) 在土壤调查样本较为稀疏的大面积高寒山区，通过把丰富的成土环境因素变量与复杂的机器学习方法相结合，模型解释度达83% 以上，能够实现较为准确的土壤养分预测制图，捕捉土壤空间变异。

2) 青海省表层土壤养分在整体上呈现东南–西北均匀递减趋势。植被能有效地指示青海省表层土壤养分全氮(TN)、全磷(TP)、全钾(TK)空间分布模式，其中年均降水量、地表温度是影响青海省表层土壤TN空间模式的重要因素；地表覆被、海拔和地表温度等环境因子对表层土壤TP的空间变异占主导作用；年均降水量、昼夜温差对表层土壤TK的空间变异占主导作用。

3) 土壤养分预测结果高等级出现在南部的玉树、果洛、黄南和东部的湟水谷地地区，低等级主要分布在柴达木盆地、可可西里、海南州中南部，全省养分等级呈中上等级以上，占全省面积的81%。由此建立的土壤养分管理分区，可为全省土壤养分的监测与农牧业区划的管理提供参考依据。