基于主成分聚类法的典型黑土区土壤地球化学分类

刘凯，戴慧敏，刘国栋，宋运红，梁帅，杨泽

(1.中国地质调查局 沈阳地质调查中心，辽宁 沈阳 110034;2.自然资源部 黑土地演化与生态效应重点实验室，辽宁 沈阳 110034;3.辽宁省黑土地演化与生态效应重点实验室，辽宁 沈阳 110034)

0 引言

东北黑土地耕地面积占全国面积的27%，粮食总产量占全国的1/4，是我国十分重要的商品粮基地，也是保障国家粮食安全的“压舱石”。2003～2020年，中国地质调查局在我国平原区持续开展了多目标区域地球化学调查工作[1]，其中在东北平原区已完成调查面积约48万km2，采集了60余万件土壤样品，获取了海量的土壤地球化学数据，为开展黑土区地球化学研究提供了庞大的数据集[2]。目前，该数据集在黑土区碳储量计算[3-5]、土地质量评价[6]、关键带物质循环[7-8]等方面得到了广泛应用，但关于黑土区常量元素的研究尚未见报道。

常量元素是土壤化学成分中最主要的组成部分，能有效反映土壤的粒度组成、矿物成分等属性[9]。同时，在土壤的形成演化过程中，常量元素也扮演着重要角色[10-12]。更重要的是，由于常量元素性质相对稳定，通常可以反映成土母质的化学组成，进而可以利用常量元素对成土母质进行有效预测[13-14]。

地质过程必然会在土壤地球化学数据中留下痕迹[15]，多种分析方法的结合更有助于揭示元素地球化学与土壤和质体间的相互关系[16]。主成分分析和聚类分析都是研究土壤地球化学的常用统计方法，主成分分析可以通过降维技术把多个变量化为少数几个主成分，这些主成分仍能够反映原始变量的大部分信息[17]；聚类分析可以将性质相似的样品归为一类，从而达到分类的目的。主成分聚类法则能将两种方法的优势进行融合，先通过主成分分析获得样品的主成分得分值，再利用得分值进行聚类分析。目前一些学者将主成分聚类法应用于土地质量评价[18]、土样识别[19]以及土壤管理分区[20]等方面，但在土壤地球化学分类中的应用仍然较少。

本文基于多目标区域地球化学获得的土壤常量元素数据，尝试利用主成分聚类法，对东北典型黑土区进行地球化学分类研究并探讨其地质意义。

1 研究区概况

目前关于东北黑土存在“广义黑土区”、“黑土区”、“典型黑土区”、“黑土地”等多种说法[21-22]，其定义和分布范围差异较大。本文研究对象主要为中国土壤发生学分类中的黑土区，但由于其空间上分布分散，与草甸土等土壤相互交错，为保证研究区的完整性，根据中国1∶100万土壤图(来源为http://westdc.westgis.ac.cn)绘制了黑土的轮廓作为本次研究区(图1)，总面积约7万km2。

图1 研究区区域位置(a)及土壤类型分布(b)Fig.1 Regional location(a) and soil typedistribution(b) of the study area

研究区位于松嫩平原东部、北部的山前台地和高平原区，北起黑龙江省嫩江县，南至辽宁省昌图县，沿哈尔滨—北安、哈尔滨—长春铁路沿线形成一条完整的黑土带。该区年降雨量在500～600 mm，绝大部分集中于7～9月，干燥度≤1。年平均气温0.5～6 ℃，土壤冻结时间为120～200 d。研究区土壤类型以黑土为主，其次为草甸土、暗棕壤，也包含有少量的黑钙土、沼泽土和火山灰土。研究区尚未编制精度较高的成土母质图，不同成土母质的分布尚不清晰，文献中多简单叙述为以第四纪更新世砂砾、黏土层和全新世砂砾、黏土层为主[23]。

2 材料和方法

2.1 数据来源

本次研究采用的土壤地球化学数据全部来源于多目标区域地球化学调查[2]。土壤样品获取采用双层网格化土壤测量方法，分别采集表层(0～20 cm)和深层(150 cm以下)土壤，采样质量为1 kg，样品在自然风干后过20目尼龙筛。表层土壤采样密度1点/km2，深层土壤采样密度1点/4 km2，每4个原始样品组合后测试分析54项元素和指标，因此每个表层测试样品代表4 km2范围，每个深层测试样品代表16 km2范围。测试元素和指标共54项，其中常量元素SiO2、Al2O3、Fe2O3、MgO、CaO、Na2O、K2O采用X荧光光谱仪测试完成，测试准确度和精密度均达到相关规范[24]要求。在剔除个别异常点后，研究区内共选取22 888个表层土壤样品和5 754个深层土壤样品用于本次研究。

2.2 数据处理方法

2.2.1 描述性统计

元素含量的平均值、中位数、标准差等利用SPSS软件进行计算。富集系数(q)代表元素在表层土壤和深层土壤中含量的比值，在ArcGIS中按空间对应关系将1个深层样品化学属性赋予4个表层土壤样品，进而计算每个表层土壤的q值。变异系数(CV)是标准差与平均值的比值，用以代表元素含量的离散程度。

2.2.2 主成分分析

主成分分析是常用的多元统计方法之一，主要用于多元数据的降维，本文使用SPSS 软件进行主成分分析操作。主成分分析的前提是参与分析的各变量间具有相关性，图2显示各常量元素之间多具有显著相关性，KMO检验系数为0.6，Bartlett’s检验结果P<0.001，说明常量元素之间存在共线性，可进行主成分提取。

图2 相关系数矩阵Fig.2 Correlation coefficient matrix diagram

2.2.3 K均值聚类

聚类分析是无监督学习的一个重要工具，包括多种具体的聚类方法，如层次聚类分析、基于模型的聚类分析和模糊聚类等[25]。K均值聚类是最典型的聚类方法之一，其基本思想是以空间中k个点为中心对样品进行聚类，将最靠近中心的样品归为一类。该算法的最大优势在于简洁和快速，自Macqueen[26]提出以来，该方法已成为大量数据分析中最流行的算法之一，对于区域地球化学调查获取的海量样品最为适宜。

K均值聚类由于要手动输入聚类数目，因此其难点之一就是确定最优分类数(K值)[27]。通常，当对数据集的属性非常了解时，可以靠经验来确定分类数，但这对数据处理者的知识储备提出了非常高的要求，往往难以实现。因此，大多数情况下人们通过计算不同分类数的某些特征，进行比较后来决定最优分类数。安光辉等[28]利用一致性检验和区间差异显著性检验来评价不同分类数目效果，但该方法的主观性较大。Rousseuw[29]于1987年提出利用轮廓系数来判断分类数，轮廓系数结合了内聚度和分离度两种因素，可以用来在相同原始数据的基础上评价不同算法或者算法不同运行方式对聚类结果所产生的影响，目前得到较广泛的应用。轮廓系数计算公式如下：

式中：ai为第i个对象到其所属簇中所有对象的平均距离；bi为该对象到所有非所属簇中对象的平均距离。轮廓系数在-1～1间变化，越趋近于1代表内聚度和分离度都相对较优，聚类效果相对较好。因此，本文利用平均轮廓系数对分类结果进行评价。K均值聚类和轮廓系数计算利用Matlab软件实现。

3 结果与讨论

3.1 常量元素含量特征

表1为研究区表层和深层土壤的常量元素统计数据，其元素含量高低趋势一致，均为SiO2>Al2O3>Fe2O3>K2O>Na2O>CaO>MgO，其中SiO2、Al2O3和Fe2O3含量占80%以上，具有明显的硅铝土特点。从元素含量分布形态看，除CaO以外的其他元素偏度在1附近，基本符合正态分布，而CaO表层和深层偏度值分别为6.25和4.30，呈右偏形态，说明存在一定数量的极大值。

表1 黑土区常量元素参数统计Table 1 Statistical table of major elements in black soil area

富集系数q代表了元素在土壤表层和深层的变化程度，当q在0.9～1.1，说明元素含量在表层和深层变化不明显，q>1.1时说明元素在表层富集，q<0.9时说明元素在表层贫化而在深部富集。图3所示，SiO2的q值全部在0.9～1.1之间，说明在垂向上最为稳定，Al2O3、K2O、MgO、Fe2O3和Na2O的表深土壤含量比分布范围略大，但仍以0.9～1.1为主，CaO则表现为在表层更富集。变异系数CV代表了元素含量的离散程度，CaO表层和深层CV值分别为0.53和0.47，为中等强度变异[30]，说明CaO含量分布较不均匀。其他常量元素CV值均小于0.3，说明其空间变异程度较低。

图3 常量元素富集系数箱线图Fig.3 Box diagram of enrichment coefficient of major elements

以上结果表明，表层和深层土壤的常量元素整体特征相似，体现了表层土壤对成土母质的继承性[31]，而表层土壤常量元素变异系数略大，则可能是成土过程中元素的迁移和再分配造成的[32]。由于地球化学分类更突出元素含量高低和空间分布的差异性[16]，因此本次选用表层土壤用于地球化学分类研究。

图4为研究区表层土壤常量元素的空间分布图，插值方法为反距离权重法。SiO2的高值区主要分布于黑土区南部和东北部，Al2O3和Fe2O3的空间分布趋势整体上与SiO2相反。K2O和Na2O分布特征类似，在黑土区南部和西北部为高值区。CaO和MgO分布特征基本一致，表现为西高东低的分布格局。不同元素的空间分布特征有一定的相似性，但差异性也很明显，因此难以直接利用单个常量元素进行地球化学区的划分，必须借助于主成分分析、聚类分析等方法来挖掘控制元素分布特征的内在因素[33]。

图4 东北黑土区常量元素含量空间分布(各元素含量区间取值为四分位数，数据单位为%)Fig.4 Spatial distribution of major elements in black soil area of Northeast China(interval value of element content is quartile and data unit is %)

3.2 主成分分析结果

本文采用累积贡献率法来判断主成分的最优数量。表2显示前3位主成分的特征值大于1，分别解释46.708%、21.056%和15.289%的数据变异，因此本研究最终选取3个主成分(PC1、PC2、PC3)，提取后的主成分共可解释83.05%的数据变异。

表2 主成分特征值及方差贡献率Table 2 Principal component characteristics and variance contributions

图5显示PC1与Al2O3、Fe2O3、MgO和CaO为正相关，指示PC1代表黏土矿物和铁铝氧化物。PC2与CaO及MgO呈正相关，表明PC2代表了碳酸盐矿物。PC3与K2O、Na2O相关性最强，说明PC3代表长石类矿物。

图5 主成分分析得分Fig.5 Principal component analysis score diagram

3.3 K均值聚类结果

利用样品的主成分得分值进行K均值聚类。首先将分类数K值分别设置为2～8，并利用Matlab中evalclusters函数计算对应的平均轮廓系数。结果显示，当分类数为5时，平均轮廓系数达到最大值0.57(图6)，因此，将研究区内数据分为5类较为合理，使用Turkey多重比较法检验各类别间差异的显著性，结果显示各类别常量元素含量具有显著性差异(P<0.05)，各类别样品数及元素平均含量见表3。

图6 不同分类数及对应的平均轮廓系数Fig.6 Mean silhouette coefficients of differentclassification numbers

表3 各类别样品常量元素平均值统计Table 3 Statistical table of mean values of major elements in various classifications %

为研究分类结果的实际内涵及应用价值，绘制了样品类别的空间分布(图7)。

图7 样品类别空间分布Fig.7 Spatial distribution map of sample categories

第Ⅰ类样品数最少，仅为18个，但其元素含量特征明显，SiO2含量最低，平均值仅56.15%，而Fe2O3、MgO、CaO含量最高(表3)，常量元素含量特征与基性岩较一致。第Ⅰ类样品空间上位于五大连池玄武岩区，土壤类型为火山灰土。

第Ⅱ类样品与第Ⅲ类样品元素含量特征较为类似，其样品数量也占绝对优势，分别占全区的28.5%和43.8%。第Ⅱ类样品的SiO2含量(平均值64.66%)略高于第Ⅲ类样品(平均值63.07%)，反映其相对含有更高的砂质成分。在空间上，第Ⅱ类样品主要分布在小兴安岭西麓丘陵地区，海拔普遍大于200 m(图1)，该地区在第四纪时期以隆升为主，地形切割较剧烈，下部的砂层和砂砾石层多出露近地表，土壤类型以暗棕壤为主。

第Ⅲ类样品数量最多，占全区的43.8%，其常量元素特征与典型黑土的成土母质——黄土状亚黏土[34]基本一致，广泛分布在松辽平原东北部的高平原区，土壤类型多为黑土和草甸土。

第Ⅳ类样品仅778件，常量元素表现为含有较丰富的CaO和MgO，其平均含量分别达5.04%和1.7%，代表土壤中碳酸盐矿物含量较高。空间上沿个别河流分布，并均分布在隆起抬升的陡坡一侧，与黑钙土的分布吻合度较高。

第Ⅴ类样品SiO2、K2O、Na2O含量较高，说明土壤中含有较多的石英和长石类矿物。在空间上，第Ⅴ类样品主要分布在松花江南部各河流两岸及一级阶地。

3.4 地球化学分类结果与地质单元的关系

土壤中元素地球化学特征受成土母质、气候、地形、生物、时间等因素综合控制，而以往研究表明，土壤常量元素特征与成土母质的关系最为密切[32,35]，土壤常量元素的空间分布与成土母质和地质单元往往具有较高的吻合度[36]。

本次聚类结果中，除第Ⅰ类样品与玄武岩区完全对应外，其他4类样品在空间上与地质单元的对应性并不明显(图7、图8)，其统计学结果也得出同样结论(图9)。原因可能包括两个方面：一方面，地质图的比例尺与土壤地球化学调查尺度不匹配，在以往开展的基础地质调查工作中，第四系广泛覆盖的平原区并不是调查重点，东北平原区目前编制的第四系地质图仅为1∶100万比例尺，精度较低；另一方面，同一地质单元的地表岩性差异显著，第四纪地质单元通常由多个沉积地层构成，如本区冲湖积物通常上部为亚黏土，下部则为砂层、砂砾层，受地表侵蚀、沉积、风化等作用影响，同一地质单元在不同地区出露岩性有较大差异，因此，地质单元与地球化学的对应关系并不完全一致[37-38]，也是造成二者对应性低的一个原因。

图8 典型黑土区第四纪地质Fig.8 Quaternary geological map of typical black soil area

图9 不同类别所包含地质单元的累积频率Fig.9 Cumulative frequency map of geologicalunits in different classifications

由于区域多目标地球化学调查为1∶25万比例尺，采样密度大，因此地球化学分类结果比地质图能更清晰地反映成土母质的空间分布情况。综合地理、地质和土壤特征，可将5个类别分别命名为：Ⅰ—五大连池玄武岩母质发育的火山灰土；Ⅱ—小兴安岭山前丘陵冲洪积物母质发育的暗棕壤；Ⅲ—松辽平原东部高平原低钙黄土状土母质发育的黑土；Ⅳ—松辽平原东部高平原边缘富钙冲湖积母质上发育的黑钙土；Ⅴ—松辽平原东南部冲积物母质上发育的黑土。该结果较客观地反映了典型黑土区的成土母质分布情况。

3.5 地球化学分类结果反映的黑土生态环境问题

地球化学分类结果表明，松花江南部土壤中砂质含量较北部明显偏高(图7)，原因一方面是由于南部成土母质以河流冲积物、冲洪积物和残坡积物为主(图8)，普遍含有较高的砂质成分，但另一方面可能也反映出局部地区的黑土存在沙化现象。近期研究表明，中国东北部地区已成为沙尘暴的多发区，尤其是松花江、嫩江流域的局部草地退化和沙化地区[39]。受季风和沙尘暴作用影响，每年都有大量细粒沙尘由西向东吹至黑土区。以哈尔滨地区为例，现代尘暴中干沉降和湿沉降粉尘的REE模式和Sr-Nd同位素组成都表明，沙尘来源于科尔沁沙地和浑善达克沙地[40-41]，这些沙尘的加入会增加黑土中砂质含量，从而造成土壤沙化。在本次研究结果中，松花江南部冲湖积母质上发育的黑土中目前也具有较高的SiO2含量，指示这些地区土壤可能已发生沙化，未来应作为黑土区的生态问题加以关注。

4 结论

1)典型黑土区常量元素SiO2、Al2O3和Fe2O3含量占80%以上，具有明显的硅铝土特点，表层土壤和深层土壤的常量元素比以0.9～1.1为主，说明土壤常量元素主要受控于成土母质。

2)利用主成分聚类法将黑土区土壤样品划分为5类最为合理，各类别地球化学特征差异显著，其空间分布能较清晰地反映出成土母质的空间分布情况。根据地理、地质及土壤特征分别命名为：I—五大连池玄武岩母质发育的火山灰土；Ⅱ—小兴安岭山前丘陵冲洪积物母质发育的暗棕壤；Ⅲ—松辽平原东部高平原低钙黄土状土母质发育的黑土；Ⅳ—松辽平原东部高平原边缘富钙冲湖积母质上发育的黑钙土；Ⅴ—松辽平原东南部冲积物母质上发育的黑土。

3)松花江南部冲湖积母质发育的黑土局部具有较高的SiO2含量，指示该地区黑土已出现沙化生态问题，在黑土地保护中应给予关注。