土壤磁化率是环境作用和成土过程的产物。磁化率的大小与有机质、土壤含水量、氧化铁含量及形态等土壤属性有关，广泛应用于环境变化历史的重建^[1]。在大空间尺度上，土壤磁化率是母质、地形、气候和植被等要素对土壤属性影响的结果，其因环境因子的空间变化而呈一定规律地变化^[2]。

目前，国内外学者对不同地区沉积物和现代土壤的磁化率与环境因素之间的相关性进行了大量报道，相关研究重点围绕土壤磁化率及相关磁学指标与气候之间相关性^[3-4]。中国表层土壤磁化率的研究表明，磁化率与年均降水量和年均气温呈显著相关性，因此可用于古气候重建工作。中国黄土高原的研究普遍认为，土壤磁化率随着降水和气温的增加而增加，因为降水和气温的增加促进土壤铁氧化物形成^[5]。但是，这与美国阿拉斯加黄土的研究结果相反 ^[6]。中国红壤区也发现，土壤磁化率随降水增加而减少的现象^[7]。

不同学者对磁化率与气候之间的相关性认识不一致，因为气候仅是土壤磁化率变化的影响因素之一。在海南岛的研究发现，玄武岩发育的土壤磁化率高于其他母岩发育的土壤1~2个数量级^[8]。在英格兰的研究发现，地形和排水是影响土壤磁化率的主要因素，滞水条件下Fe³⁺还原为Fe²⁺导致土壤磁性减弱。还有研究发现，植被通过向土壤输送有机质影响土壤磁化率，富含有机质的土壤铁磁性矿物溶解，高矫顽力矿物相对富集^[9]。因此，不同区域的土壤发育和磁化率的主控因素可能存在差异，磁化率是多种环境要素综合作用的结果。

青藏高原是世界上最年轻且水平地带性和垂直地带性紧密结合的自然地理单元，影响区域与全球的能量和水分循环；因经历6次大地质构造活动，母岩类型多样，是全球变化对地表影响最显著的地区之一。本文系统调查了青藏高原范围内254个表层（发生层A层）土壤样品，旨在理解该区土壤磁化率空间变化格局，以及多种环境因素对表层土壤磁化率的可能影响机制，以期土壤磁化率能作为其他环境变量的代用指标。假设：（1）土壤磁化率的空间变异与气候、土壤等环境因素呈显著相关性；（2）植被指数是气候、地形和土壤等要素共同作用的结果，能更好地反映磁化率。

青藏高原分布范围大致为24.8—39.2°N，74.9 —104.3°E，海拔864~5 712 m，属于高寒山地气候，是全球气候变化最敏感的地区之一。采样点的年平均温度分布在–4~12℃、年平均降水量分布在40~790 mm，南部边缘为古热带植物区，东南部为中国-喜马拉雅森林植物亚区，柴达木盆地为亚洲荒漠植物亚区，腹地为青藏高原植物亚区^[10]，坡度为0~25°，主要土壤类型为新成土和雏形土^[11]。

根据典型土壤类型面积百分比估测主要土壤类型样点数，尽量选择不易积水、不易受侵蚀和堆积的地形且远离人为干扰的样点，采集发生层A层的土壤。共采集254个表土样品，这些样本能够代表研究区主要的土壤类型。研究区中北部因道路不可达，故没有布设样点。

图 1 研究区表层采样点分布 Fig. 1 Distribution of sampling points in the surface soil layer of the research area

采集的土样在室内经自然风干、去杂和研磨过不同孔径筛后，用于土壤理化属性的测定。其中，颗粒组成采用吸管法测定；pH采用电位法测定；有机碳采用重铬酸钾外加热法测定；全磷采用钼锑抗比色法测定；阳离子交换量（CEC）采用乙酸铵（pH7.0）-EDTA交换法测定^[12]。本文使用P-XRF（Amptek，Bedford，MA，USA）测定了土壤样品中K、Ca、Fe、Si和Ti等主要元素浓度，测量时间为2 min。土壤游离氧化铁和pH分别采用傅立叶变换红外光谱仪（Nicolet iS50，Thermo Scientific，USA）和pH计（METTLER TOLEDO，USA）测定。

在空气干燥的条件下，将土壤样品放置在2 cm的正方体塑料盒中，进行磁性测定。使用卡帕桥多频磁化率仪（AGICO MFK1 FA）在976 Hz（低频磁化率，χ_lf）和15 616 Hz（高频磁化率χ_hf）的双重频率下测量特定质量（6 g）的磁感应强度（magnetic susceptibility，χ）。卡帕桥系列是测量磁化率的主要仪器，测量精度要较Bartington磁化率仪高，预计会成为下一代磁化率仪器的主流^[13]。χ_lf为低频磁化率，χ_hf为高频磁化率，χ_fd为频率磁化率，χ_fd%为百分频率磁化率。χ_fd = χ_lf – χ_hf，χ_fd% =（χ_lf – χ_hf）/χ_lf × 100。本研究随机选取7个样品，分别用AGICO MFK1 FA（低频磁化率χ_lf-K，高频磁化率χ_fd%_-K）和Bartington（低频磁化率χ_lf-B，高频磁化率χ_fd%_-B）磁化率仪器测定磁化率，发现二者测得磁化率数值之间的关系为：χ_lf-B = 1.2χ_lf-K + 0.02，R² = 0.999；χ_fd%_-B=10^–3χ_fd%_-K + 10^–4，R²=0.998。所有的磁化率测量均在中国科学院南京地理与湖泊研究所磁化率分析实验室进行。

根据相关研究结论，本研究收集了11种驱动土壤演化的环境变量（表 1）。气候变量为1970—2020年的均值，从WorldClim2数据库中获得^[14]，分辨率为1 km。地形变量来自90 m SRTM DEM（Consortium for Spatial Information（CGIAR-CSI），2018），在SAGA GIS^[15-16]中计算得出。表层地质（土壤母质）用30 m水平分辨率的Landsat8 ETM + Band7（2.08~2.35um的短波红外）和黏土矿物指数（Band5/Band7）表示^[17]。植被使用2000—2017年每月最大归一化植被指数（NDVI）的平均值、最大值和最小值来表示。

表 1 Table 1

表 1 青藏高原气候、母质、地形和植被等环境协变量的分辨率和来源 Table 1 Scale/resolution and source for environmental covariates of climate, parent material, terrain and vegetation 类别Category 协变量Covariate 描述Description 分辨率Scale/resolution 来源Source 气候Climate MAT 年均气温 1 km RESDC① MAP 年均降水量 1 km RESDC Max MAP 最高年均降水量 1 km RESDC Min MAP 最低年均降水量 1 km RESDC Max MAT 最高年均气温 1 km RESDC Min MAT 最低年均气温 1 km RESDC Aridity 干旱系数 1 km RESDC 母质Parent material Rock type 岩石类型 1︰2.5 M RESDC 地形Terrain DEM 海拔 90 m Jarvis et al. [18] Slope 坡度 90 m Conrad et al. [15] Slopelength 坡长 90 m Conrad et al.[15] Curvature 曲率 90 m Conrad et al.[15] 植被Vegetation NDVImean 归一化植被指数平均值 30 m USGS② NDVImax 归一化植被指数最大值 30 m USGS NDVImin 归一化植被指数最小值 30 m USGS ①RESDC：中国科学院资源环境科学与数据中心（Data Center for Resources and Environmental Sciences Chinese Academy of Sciences）； ②USGS：美国地质调查局，https://www.usgs.gov/

表 1 青藏高原气候、母质、地形和植被等环境协变量的分辨率和来源 Table 1 Scale/resolution and source for environmental covariates of climate, parent material, terrain and vegetation

采用Pearson方法分析了254个样品的磁化率与环境变量以及土壤属性的相关性。使用对数正态变换将χ_lf数据转换为卡方分布P＜0.05的正态分布数据log₁₀χ_lf。用已测量的土壤磁化率与环境变量构建随机森林（RF）模型，预测χ_lf和χ_fd%的空间分布图。RF的预测性能通过平均误差（MAE）、均方根误差（RMSE）和决定系数（R²）的10折交叉验证来评估。使用R软件（版本3.5.1，http://cran.r-project.org/）的“randomForest”包实现土壤磁化率制图。采用主成分分析和逐步回归的方法分析土壤磁化率的主控因子。采用Excel 2017进行数据整理，R 3.6.3进行基本特征分析和绘图。

青藏高原表层土壤磁化率变幅较大，χ_lf变化范围为3.0 × 10^–8~5.979 × 10^–6 m³·kg^–1，均值为6.23 × 10^–7 m³·kg^–1；χ_fd%变化范围为0~17.3%，均值为5.2%（表 2）。χ_lf的极大值位于拉萨市。除了结晶盐风化物以外，不同母质类型之间，表层土壤χ_lf无显著差异；χ_fd%差异显著：黄土和砂页岩风化物 > 风沙、冰碛物和结晶盐风化物，其他母质类型之间无显著差异（图 2）。不同土地利用类型之间，表层土壤χ_lf无显著差异，χ_fd%差异显著：森林和旱地 > 草地 > 荒地（图 3）。

表 2 Table 2

表 2 本文数据集的统计信息 Table 2 Summary statistics for the datasets discussed in the main text 变量 Parameters 均值 Mean 误差 S.D. 中位数 Median 最小值 Min. 最大值 Max. 偏离度 Skew 峰度Kurtosis 磁化率χlf/（10–8 m3·kg–1） 62.3 64.7 42.0 3.0 597.9 3.9 23.8 百分频率磁化率χfd%/% 5.2 4.2 4.4 0.0 17.2 0.1 –0.5 归一化植被指数平均值NDVImean 0.5 0.2 0.6 0.0 0.8 –0.6 –0.9 年均降水MAP/（mm） 445.4 199.4 463.0 28.7 910.8 –0.4 –0.6 年均气温MAT/（℃） 3.6 4.3 3.0 –5.2 23.7 0.8 1.4 酸碱度pH 7.8 1.2 8.1 3.4 10.0 –1.1 0.8 有机碳SOC/（g·kg–1） 28.9 25.5 20.8 0.8 113.3 1.1 0.5 黏粒Clay/% 5.7 8.4 1.3 0.2 48.3 0.2 0.3 粉粒Silt/% 16.5 21.6 4.9 0.3 68.6 0.1 –0.0 砂粒Sand/% 12.9 15.3 6.3 0.9 90.2 0.2 0.5 铁Fe/（g·kg–1） 57.0 12.0 57.6 28.7 86.1 –0.1 0.1 游离氧化铁Free iron oxide/（g·kg–1） 22.8 8.1 22.5 8.7 48.9 0.8 0.9

表 2 本文数据集的统计信息 Table 2 Summary statistics for the datasets discussed in the main text

注：方框和竖线代表每个箱子中的最小值、第25百分位、中值、第75百分位和最大值。框中的黑色实线表示平均值。黑色实点为异常值，表示为位于箱线图边缘之外的数据点。竖线下面的数字为样本大小。竖线上方不同小写字母表示不同母质在P < 0.05水平上存在显著差异。下同。  Note: Boxes and whiskers represent the minimum, 25th percentile, median, 75th percentile, and maximum values in each bin. The black solid lines in the boxes represent the mean values. Black solid points are outliers that are defined as data points located outside the whiskers of the boxplots. The numbers below the whiskers are sample sizes. Different lowercase letters above the whiskers denote significant difference in different parent materials at the P < 0.05 level.The same as below. 图 2 对比不同母质之间（a）χlf和（b）χfd%的差异 Fig. 2 Comparison of(a)χlf and(b)χfd% in different parent materials

图 3 不同土地利用类型之间χlf（a）和χfd%（b）的差异 Fig. 3 Comparison ofχlf(a)andχfd%(b)in different land uses

基于随机森林（RF），χ_lf的预测精度RMSE和R²分别为51.7和0.13。χ_fd%的预测精度RMSE和R²分别为3.4和0.30。χ_lf和χ_fd%的半变异函数模型最好通过高斯模型拟合（图 4），其中参数通过限制最大似然估计。χ_lf和χ_fd%的块金值与基台值比率大于50%（图 4），表明χ_lf和χ_fd%的空间依赖性结构较弱，56%~62%的χ_lf和χ_fd%变异性由随机因素组成。RF在χ_lf（R² = 0.13）和χ_fd%（R² = 0.30）预测方面略优于地理统计空间插值方法（χ_lf的R² = 0.09和χ_fd%的R² = 0.30），这主要归因于弱空间依赖结构。因此，本文用已测量的土壤磁化率（χ_lf和χ_fd%）与环境变量之间的关系，构建RF模型，用于绘制χ_lf和χ_fd%的空间分布图（图 4）。

图 4 基于随机森林模型模拟的表层土壤（发生层A层）χlf（a）和χfd%（b）空间分布图 Fig. 4 Maps of χlf(a)and χfd% (b)predicted by the random forest analyses in topsoils(occurrence layer A)

宏观上，东南部的χ_lf和χ_fd%均显示较高的值（图 4）。χ_lf空间分布呈由东南向西北降低的趋势，高值区位于东南部边缘，低值区位于腹地。χ_fd%空间分布规律与χ_lf相似，高值区位于东部偏南，低值区位于西部。小尺度地形上，χ_lf和χ_fd%的空间分布显示许多细节，主要是由于环境变量的变化对RF模型预测的影响，比如在河谷地带χ_lf和χ_fd%的值较周边海拔较高的地点高。

青藏高原MAP、MAT、DEM和NDVI具有明显的空间分布规律（图 5）。MAP呈弧形带状分布，由东向西递减；MAT大尺度地形区呈由南往北递减，小尺度地形区呈与DEM相匹配的分布规律——随海拔升高而递减。NDVI大尺度地形区上呈由东南向西北递减趋势，小尺度地形区呈与DEM相匹配的分布规律，南部边缘为最高，东南部其次，腹地最低，与χ_lf和χ_fd%的空间分布规律吻合。

图 5 青藏高原各环境因子空间分布 Fig. 5 Distribution of environmental factors in the Qinghai-Tibet Plateau

通过对土壤磁化率与对应环境因子和土壤属性进行Pearson相关分析（图 6）发现，χ_lf对数转化后的Log₁₀χ_lf与土壤母质代表性物质（Fe和Ti）的相关性最高，与SOC和TN呈显著负相关；与小尺度地形区上的坡度（Slope）和坡长（Slopelength）呈显著正相关；与气候（MAP和MAT）无显著相关。除了与土壤属性和地形显著相关外，χ_fd%还与气候（MAP和MAT）和植被指数（NDVI）显著相关，相关性土壤属性和植被指数 > 气候因子 > 地形因子。

图 6 土壤磁化率与环境影响因素之间的相关性 Fig. 6 Correlation between soil magnetic susceptibility and environmental factors

土壤磁化率与多种环境因子的相关系数矩阵（图 6）在反映土壤磁化率与各因子之间的关系的同时，也说明多种因子之间存在多重共线性关系。其中，多个环境因子和土壤属性之间相关系数的绝对值大于0.5，且均为极显著水平。为了筛选出土壤磁化率的主控因子和确定各因子对土壤磁化率的贡献，采用主成分分析和逐步回归的方法分析土壤磁化率的主控因子。

主成分分析结果显示，前5个主成分保留了大部分原影响因子的信息，其特征值> 1，PC1、PC2和PC3的贡献率均大于10%，PC4的贡献率9%，PC5的贡献率6%（表 3）。第一主成分（PC1）主要受植被归一化指数的影响，命名其为植被主成分；第二主成分主要受土壤颗粒（Clay、Silt和Sand）的影响，命名其为土壤物理属性主成分；第三主成分主要受DEM的影响，命名其为大尺度地形区主成分；第四主成分主要受坡度和坡长的影响，命名其为小尺度地形区主成分；第五主成分主要受土壤铁（Fe）和钛（Ti）的影响，二者在土壤弱风化地区一定程度上反应母质磁性强弱，命名其为母质主成分。前5个主成分的累积贡献率76%。通过表 3主成分分析结果，各主成分对土壤磁化率影响：植被主成分 > 土壤物理属性主成分 > 大尺度地形区主成分 > 小尺度地形区主成分 > 母质主成分。

表 3 Table 3

表 3 土壤磁化率环境影响因素主成分（PC）分析 Table 3 Principal component(PC)analysis of environmental factors affecting χlf and χfd% 影响因子Influence factor 主成分1 PC1 主成分2 PC2 主成分3 PC3 主成分4 PC4 主成分5 PC5 铁Fe 0.24 0.30 0.05 0.24 0.48 钛Ti 0.27 0.21 –0.09 0.23 0.49 酸碱度pH –0.26 0.18 0.26 –0.18 –0.04 黏粒Clay 0.18 0.42 –0.21 –0.16 –0.22 粉粒Silt 0.25 0.40 –0.04 –0.23 –0.18 砂粒Sand –0.25 –0.44 0.08 0.23 0.21 有机碳SOC 0.26 –0.31 –0.32 0.02 –0.08 全氮TN 0.27 –0.27 –0.33 0.03 –0.14 年均气温MAT 0.20 –0.07 0.30 –0.28 0.12 年均降水MAP 0.21 –0.17 0.23 –0.38 –0.10 海拔DEM –0.08 0.05 –0.52 0.17 –0.12 坡度曲率Curvature 0.15 0.02 0.35 0.23 –0.19 坡度Slope 0.19 0.07 0.25 0.48 –0.12 坡长Slopelength –0.12 –0.03 –0.19 –0.42 0.53 归一化植被指数均值NDVImean 0.34 –0.16 –0.01 –0.08 0.01 归一化植被指数最大值NDVImax 0.35 –0.17 0.01 –0.05 0.06 归一化植被指数最小值NDVImin 0.30 –0.21 0.16 –0.08 0.05 统计参数Statistical parameter 特征值Eigenvalues 2.49 1.50 1.35 1.23 1.04 贡献率Contribution rate/% 36 13 11 9 6 累计贡献率Cumulative contribution rate/% 36 50 61 70 76

表 3 土壤磁化率环境影响因素主成分（PC）分析 Table 3 Principal component(PC)analysis of environmental factors affecting χlf and χfd%

在主成分分析的基础上，通过逐步回归分析（引入变量P值为0.05，剔除变量的P值为0.1）筛选出更准的土壤磁化率主要调控因子。结果表明，对土壤χ_lf贡献最大的因子为Fe、SOC、MAP、DEM和Slope；对土壤χ_fd%贡献最大的因子为Fe、Ti、Silt、Sand、TN、MAP、NDVImean和NDVImax。

在区域尺度下，土壤磁化率主要受母质、气候、地形、植被和土壤的影响^[19]。这与本研究中土壤磁化率受多种环境因子影响的结果一致。本研究结果表明，土壤磁化率与植被、土壤物理属性、地形和母质之间存在极显著相关性。

土壤磁性矿物的来源首先是母岩的原生矿物，其次是成土过程中产生的次生矿物，前者一般为砂粒或粗粉粒的粒度级；后者一般为稳定单畴（Stable single domain，SSD）粒度或更小，且有黏土的粒度级^[20]。浙江红壤χ_lf随粒径减小而减小，说明原生矿物残留对浙江红壤磁化率的贡献为主，次生磁矿物为辅^[20-21]。相反，云贵高原土壤磁化率χ_lf随粒径减小而增大，说明次生矿物对其磁化率变化起主要贡献^[20]。本研究中，χ_lf和χ_fd%均随粒径减少而增大（图 6），说明次生矿物对其磁化率变化同样起主要贡献。

国内外学者已经开展大量表土磁化率与现代气候因子关系的研究^[22]。本研究中，χ_fd% 与年均降水量和年均气温显著正相关，与黄土高原及周边地区的磁化率与气候的相关性一致，与长江以南地区的磁化率与气候的相关性相反^[23]。已有研究表明，年均降水量低于1 000~1 500 mm时，随年均降水量和年均气温的增加，土壤铁氧化物增加，从而促进土壤磁化率增加；而中国南方表层土壤磁化率随年均降水和年均气温的降低而增高，重要原因之一可能是高温多雨的气候使得土壤水分饱和，土壤中的磁性矿物易被还原和分解，同时，湿热环境中的强磁性磁铁矿易于转化为弱磁性赤铁矿^[23-25]。χ_lf与气候的相关性不显著，可能青藏高原次生磁化率对青藏高原的土壤磁性变化起主要贡献（图 6），或者采样点的实际气候受到海拔的影响，亦或者青藏高原的气象站点稀疏，导致气象资料误差。

在英格兰的研究发现，局部尺度上，地形、土地利用和植被影响土壤形成和铁磁性矿物的累积^[26-30]。本研究中，海拔调控热水在垂直方向上的再分配影响植被^[29]，进而间接影响土壤磁化率；土壤磁化率随坡度增加而增加，坡度越大有利于排水，与已有研究一致。该研究表明，地形主要通过调控土壤排水和水热条件再分配，进而影响土壤磁化率^[27]。

自然土壤的表土层，一般具有较高的磁性，与生物成磁作用有关^[22]。本研究中，χ_lf均值和χ_fd%均值为森林 > 草地 > 荒地，与已有大多数研究结果一致，且与NDVI值变化趋势相吻合（图7）。土壤有机碳为铁还原菌提供了能量，这些细菌将Fe³⁺还原为Fe²⁺，过量的Fe²⁺与铁氢化物反应生成SP或SD磁铁矿^[27]。此外，本研究中旱地与森林无显著差异，说明本研究区人为活动使表土磁化率增强，与澳大利亚东南部的研究结果一致^[2]。人类活动强烈区域中大多数空气中的氧化铁均以聚集的磁铁矿纳米颗粒的形式存在^[28]，这些磁性颗粒随着雨水或者重力作用沉降到表土，增强表土磁性。

吕厚远等^[23]研究认为，土壤磁化率存在区域尺度上的分布规律。本研究表明，青藏高原土壤χ_lf空间分布呈东南向西北降低的趋势，高值区位于东南部边缘，低值区位于中部。χ_fd%空间分布规律的趋势与χ_lf相似，高值区位于东部偏南，低值区位于西部。这主要是由于区域尺度上的土壤磁化率受多种环境因子的影响，本研究中环境因子空间分布趋势（图 6）以及各要素之间的相关分析（图7）均证明这一点。

值得注意的是，本研究χ_lf和χ_fd%空间分布规律与青藏高原植被类型分区吻合，χ_lf和χ_fd%值在南部边缘最高，东南部和柴达木盆地其次，腹地最低，分别对应南部边缘古热带植物区，东南部中国-喜马拉雅森林植物亚区，柴达木盆地亚洲荒漠植物亚区，腹地青藏高原植物亚区^[10]。植被对土壤次生磁性矿物的形成产生重要影响（图 3b），比如，森林χ_fd%显著大于荒地，因为森林土壤的透水性和持水性较适中，而荒地处于干旱或半干旱状态，土壤湿度对土壤矿物转化及生成有重要影响，当土壤湿度适中时，成壤过程可生成更多磁铁矿与磁赤铁矿等亚铁磁性矿物^[22]。本研究主成分分析表明，植被是影响土壤属性变异的主要因子，研究区东南部边缘为NDVI的高值区，同时是土壤磁化率的高值区，反映土壤磁化率的分布格局与NDVI关系紧密，说明磁化率能更好地反映植被空间分布。

区域尺度上，青藏高原土壤磁化率的空间分布规律为沿东南向西北方向递减的变化趋势，主要受植被、土壤属性、地形和母质的影响。青藏高原表层土壤磁化率受气候影响较小，受植被的影响更大。未来应在土壤磁化率是否能指示植被的研究上继续探索。