基于主成分聚類法的典型黑土區(qū)土壤地球化學(xué)分類

劉凱，戴慧敏，劉國(guó)棟，宋運(yùn)紅，梁帥，楊澤

(1.中國(guó)地質(zhì)調(diào)查局 沈陽(yáng)地質(zhì)調(diào)查中心，遼寧 沈陽(yáng) 110034;2.自然資源部 黑土地演化與生態(tài)效應(yīng)重點(diǎn)實(shí)驗(yàn)室，遼寧 沈陽(yáng) 110034;3.遼寧省黑土地演化與生態(tài)效應(yīng)重點(diǎn)實(shí)驗(yàn)室，遼寧 沈陽(yáng) 110034)

0 引言

東北黑土地耕地面積占全國(guó)面積的27%，糧食總產(chǎn)量占全國(guó)的1/4，是我國(guó)十分重要的商品糧基地，也是保障國(guó)家糧食安全的“壓艙石”。2003～2020年，中國(guó)地質(zhì)調(diào)查局在我國(guó)平原區(qū)持續(xù)開(kāi)展了多目標(biāo)區(qū)域地球化學(xué)調(diào)查工作[1]，其中在東北平原區(qū)已完成調(diào)查面積約48萬(wàn)km2，采集了60余萬(wàn)件土壤樣品，獲取了海量的土壤地球化學(xué)數(shù)據(jù)，為開(kāi)展黑土區(qū)地球化學(xué)研究提供了龐大的數(shù)據(jù)集[2]。目前，該數(shù)據(jù)集在黑土區(qū)碳儲(chǔ)量計(jì)算[3-5]、土地質(zhì)量評(píng)價(jià)[6]、關(guān)鍵帶物質(zhì)循環(huán)[7-8]等方面得到了廣泛應(yīng)用，但關(guān)于黑土區(qū)常量元素的研究尚未見(jiàn)報(bào)道。

常量元素是土壤化學(xué)成分中最主要的組成部分，能有效反映土壤的粒度組成、礦物成分等屬性[9]。同時(shí)，在土壤的形成演化過(guò)程中，常量元素也扮演著重要角色[10-12]。更重要的是，由于常量元素性質(zhì)相對(duì)穩(wěn)定，通?？梢苑从吵赏聊纲|(zhì)的化學(xué)組成，進(jìn)而可以利用常量元素對(duì)成土母質(zhì)進(jìn)行有效預(yù)測(cè)[13-14]。

地質(zhì)過(guò)程必然會(huì)在土壤地球化學(xué)數(shù)據(jù)中留下痕跡[15]，多種分析方法的結(jié)合更有助于揭示元素地球化學(xué)與土壤和質(zhì)體間的相互關(guān)系[16]。主成分分析和聚類分析都是研究土壤地球化學(xué)的常用統(tǒng)計(jì)方法，主成分分析可以通過(guò)降維技術(shù)把多個(gè)變量化為少數(shù)幾個(gè)主成分，這些主成分仍能夠反映原始變量的大部分信息[17]；聚類分析可以將性質(zhì)相似的樣品歸為一類，從而達(dá)到分類的目的。主成分聚類法則能將兩種方法的優(yōu)勢(shì)進(jìn)行融合，先通過(guò)主成分分析獲得樣品的主成分得分值，再利用得分值進(jìn)行聚類分析。目前一些學(xué)者將主成分聚類法應(yīng)用于土地質(zhì)量評(píng)價(jià)[18]、土樣識(shí)別[19]以及土壤管理分區(qū)[20]等方面，但在土壤地球化學(xué)分類中的應(yīng)用仍然較少。

本文基于多目標(biāo)區(qū)域地球化學(xué)獲得的土壤常量元素?cái)?shù)據(jù)，嘗試?yán)弥鞒煞志垲惙?，?duì)東北典型黑土區(qū)進(jìn)行地球化學(xué)分類研究并探討其地質(zhì)意義。

1 研究區(qū)概況

目前關(guān)于東北黑土存在“廣義黑土區(qū)”、“黑土區(qū)”、“典型黑土區(qū)”、“黑土地”等多種說(shuō)法[21-22]，其定義和分布范圍差異較大。本文研究對(duì)象主要為中國(guó)土壤發(fā)生學(xué)分類中的黑土區(qū)，但由于其空間上分布分散，與草甸土等土壤相互交錯(cuò)，為保證研究區(qū)的完整性，根據(jù)中國(guó)1∶100萬(wàn)土壤圖(來(lái)源為http://westdc.westgis.ac.cn)繪制了黑土的輪廓作為本次研究區(qū)(圖1)，總面積約7萬(wàn)km2。

圖1 研究區(qū)區(qū)域位置(a)及土壤類型分布(b)Fig.1 Regional location(a) and soil typedistribution(b) of the study area

研究區(qū)位于松嫩平原東部、北部的山前臺(tái)地和高平原區(qū)，北起黑龍江省嫩江縣，南至遼寧省昌圖縣，沿哈爾濱—北安、哈爾濱—長(zhǎng)春鐵路沿線形成一條完整的黑土帶。該區(qū)年降雨量在500～600 mm，絕大部分集中于7～9月，干燥度≤1。年平均氣溫0.5～6 ℃，土壤凍結(jié)時(shí)間為120～200 d。研究區(qū)土壤類型以黑土為主，其次為草甸土、暗棕壤，也包含有少量的黑鈣土、沼澤土和火山灰土。研究區(qū)尚未編制精度較高的成土母質(zhì)圖，不同成土母質(zhì)的分布尚不清晰，文獻(xiàn)中多簡(jiǎn)單敘述為以第四紀(jì)更新世砂礫、黏土層和全新世砂礫、黏土層為主[23]。

2 材料和方法

2.1 數(shù)據(jù)來(lái)源

本次研究采用的土壤地球化學(xué)數(shù)據(jù)全部來(lái)源于多目標(biāo)區(qū)域地球化學(xué)調(diào)查[2]。土壤樣品獲取采用雙層網(wǎng)格化土壤測(cè)量方法，分別采集表層(0～20 cm)和深層(150 cm以下)土壤，采樣質(zhì)量為1 kg，樣品在自然風(fēng)干后過(guò)20目尼龍篩。表層土壤采樣密度1點(diǎn)/km2，深層土壤采樣密度1點(diǎn)/4 km2，每4個(gè)原始樣品組合后測(cè)試分析54項(xiàng)元素和指標(biāo)，因此每個(gè)表層測(cè)試樣品代表4 km2范圍，每個(gè)深層測(cè)試樣品代表16 km2范圍。測(cè)試元素和指標(biāo)共54項(xiàng)，其中常量元素SiO2、Al2O3、Fe2O3、MgO、CaO、Na2O、K2O采用X熒光光譜儀測(cè)試完成，測(cè)試準(zhǔn)確度和精密度均達(dá)到相關(guān)規(guī)范[24]要求。在剔除個(gè)別異常點(diǎn)后，研究區(qū)內(nèi)共選取22 888個(gè)表層土壤樣品和5 754個(gè)深層土壤樣品用于本次研究。

2.2 數(shù)據(jù)處理方法

2.2.1 描述性統(tǒng)計(jì)

元素含量的平均值、中位數(shù)、標(biāo)準(zhǔn)差等利用SPSS軟件進(jìn)行計(jì)算。富集系數(shù)(q)代表元素在表層土壤和深層土壤中含量的比值，在ArcGIS中按空間對(duì)應(yīng)關(guān)系將1個(gè)深層樣品化學(xué)屬性賦予4個(gè)表層土壤樣品，進(jìn)而計(jì)算每個(gè)表層土壤的q值。變異系數(shù)(CV)是標(biāo)準(zhǔn)差與平均值的比值，用以代表元素含量的離散程度。

2.2.2 主成分分析

主成分分析是常用的多元統(tǒng)計(jì)方法之一，主要用于多元數(shù)據(jù)的降維，本文使用SPSS 軟件進(jìn)行主成分分析操作。主成分分析的前提是參與分析的各變量間具有相關(guān)性，圖2顯示各常量元素之間多具有顯著相關(guān)性，KMO檢驗(yàn)系數(shù)為0.6，Bartlett’s檢驗(yàn)結(jié)果P<0.001，說(shuō)明常量元素之間存在共線性，可進(jìn)行主成分提取。

圖2 相關(guān)系數(shù)矩陣Fig.2 Correlation coefficient matrix diagram

2.2.3 K均值聚類

聚類分析是無(wú)監(jiān)督學(xué)習(xí)的一個(gè)重要工具，包括多種具體的聚類方法，如層次聚類分析、基于模型的聚類分析和模糊聚類等[25]。K均值聚類是最典型的聚類方法之一，其基本思想是以空間中k個(gè)點(diǎn)為中心對(duì)樣品進(jìn)行聚類，將最靠近中心的樣品歸為一類。該算法的最大優(yōu)勢(shì)在于簡(jiǎn)潔和快速，自Macqueen[26]提出以來(lái)，該方法已成為大量數(shù)據(jù)分析中最流行的算法之一，對(duì)于區(qū)域地球化學(xué)調(diào)查獲取的海量樣品最為適宜。

K均值聚類由于要手動(dòng)輸入聚類數(shù)目，因此其難點(diǎn)之一就是確定最優(yōu)分類數(shù)(K值)[27]。通常，當(dāng)對(duì)數(shù)據(jù)集的屬性非常了解時(shí)，可以靠經(jīng)驗(yàn)來(lái)確定分類數(shù)，但這對(duì)數(shù)據(jù)處理者的知識(shí)儲(chǔ)備提出了非常高的要求，往往難以實(shí)現(xiàn)。因此，大多數(shù)情況下人們通過(guò)計(jì)算不同分類數(shù)的某些特征，進(jìn)行比較后來(lái)決定最優(yōu)分類數(shù)。安光輝等[28]利用一致性檢驗(yàn)和區(qū)間差異顯著性檢驗(yàn)來(lái)評(píng)價(jià)不同分類數(shù)目效果，但該方法的主觀性較大。Rousseuw[29]于1987年提出利用輪廓系數(shù)來(lái)判斷分類數(shù)，輪廓系數(shù)結(jié)合了內(nèi)聚度和分離度兩種因素，可以用來(lái)在相同原始數(shù)據(jù)的基礎(chǔ)上評(píng)價(jià)不同算法或者算法不同運(yùn)行方式對(duì)聚類結(jié)果所產(chǎn)生的影響，目前得到較廣泛的應(yīng)用。輪廓系數(shù)計(jì)算公式如下：

式中：ai為第i個(gè)對(duì)象到其所屬簇中所有對(duì)象的平均距離；bi為該對(duì)象到所有非所屬簇中對(duì)象的平均距離。輪廓系數(shù)在-1～1間變化，越趨近于1代表內(nèi)聚度和分離度都相對(duì)較優(yōu)，聚類效果相對(duì)較好。因此，本文利用平均輪廓系數(shù)對(duì)分類結(jié)果進(jìn)行評(píng)價(jià)。K均值聚類和輪廓系數(shù)計(jì)算利用Matlab軟件實(shí)現(xiàn)。

3 結(jié)果與討論

3.1 常量元素含量特征

表1為研究區(qū)表層和深層土壤的常量元素統(tǒng)計(jì)數(shù)據(jù)，其元素含量高低趨勢(shì)一致，均為SiO2>Al2O3>Fe2O3>K2O>Na2O>CaO>MgO，其中SiO2、Al2O3和Fe2O3含量占80%以上，具有明顯的硅鋁土特點(diǎn)。從元素含量分布形態(tài)看，除CaO以外的其他元素偏度在1附近，基本符合正態(tài)分布，而CaO表層和深層偏度值分別為6.25和4.30，呈右偏形態(tài)，說(shuō)明存在一定數(shù)量的極大值。

表1 黑土區(qū)常量元素參數(shù)統(tǒng)計(jì)Table 1 Statistical table of major elements in black soil area

富集系數(shù)q代表了元素在土壤表層和深層的變化程度，當(dāng)q在0.9～1.1，說(shuō)明元素含量在表層和深層變化不明顯，q>1.1時(shí)說(shuō)明元素在表層富集，q<0.9時(shí)說(shuō)明元素在表層貧化而在深部富集。圖3所示，SiO2的q值全部在0.9～1.1之間，說(shuō)明在垂向上最為穩(wěn)定，Al2O3、K2O、MgO、Fe2O3和Na2O的表深土壤含量比分布范圍略大，但仍以0.9～1.1為主，CaO則表現(xiàn)為在表層更富集。變異系數(shù)CV代表了元素含量的離散程度，CaO表層和深層CV值分別為0.53和0.47，為中等強(qiáng)度變異[30]，說(shuō)明CaO含量分布較不均勻。其他常量元素CV值均小于0.3，說(shuō)明其空間變異程度較低。

圖3 常量元素富集系數(shù)箱線圖Fig.3 Box diagram of enrichment coefficient of major elements

以上結(jié)果表明，表層和深層土壤的常量元素整體特征相似，體現(xiàn)了表層土壤對(duì)成土母質(zhì)的繼承性[31]，而表層土壤常量元素變異系數(shù)略大，則可能是成土過(guò)程中元素的遷移和再分配造成的[32]。由于地球化學(xué)分類更突出元素含量高低和空間分布的差異性[16]，因此本次選用表層土壤用于地球化學(xué)分類研究。

圖4為研究區(qū)表層土壤常量元素的空間分布圖，插值方法為反距離權(quán)重法。SiO2的高值區(qū)主要分布于黑土區(qū)南部和東北部，Al2O3和Fe2O3的空間分布趨勢(shì)整體上與SiO2相反。K2O和Na2O分布特征類似，在黑土區(qū)南部和西北部為高值區(qū)。CaO和MgO分布特征基本一致，表現(xiàn)為西高東低的分布格局。不同元素的空間分布特征有一定的相似性，但差異性也很明顯，因此難以直接利用單個(gè)常量元素進(jìn)行地球化學(xué)區(qū)的劃分，必須借助于主成分分析、聚類分析等方法來(lái)挖掘控制元素分布特征的內(nèi)在因素[33]。

圖4 東北黑土區(qū)常量元素含量空間分布(各元素含量區(qū)間取值為四分位數(shù)，數(shù)據(jù)單位為%)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 主成分分析結(jié)果

本文采用累積貢獻(xiàn)率法來(lái)判斷主成分的最優(yōu)數(shù)量。表2顯示前3位主成分的特征值大于1，分別解釋46.708%、21.056%和15.289%的數(shù)據(jù)變異，因此本研究最終選取3個(gè)主成分(PC1、PC2、PC3)，提取后的主成分共可解釋83.05%的數(shù)據(jù)變異。

表2 主成分特征值及方差貢獻(xiàn)率Table 2 Principal component characteristics and variance contributions

圖5顯示PC1與Al2O3、Fe2O3、MgO和CaO為正相關(guān)，指示PC1代表黏土礦物和鐵鋁氧化物。PC2與CaO及MgO呈正相關(guān)，表明PC2代表了碳酸鹽礦物。PC3與K2O、Na2O相關(guān)性最強(qiáng)，說(shuō)明PC3代表長(zhǎng)石類礦物。

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

3.3 K均值聚類結(jié)果

利用樣品的主成分得分值進(jìn)行K均值聚類。首先將分類數(shù)K值分別設(shè)置為2～8，并利用Matlab中evalclusters函數(shù)計(jì)算對(duì)應(yīng)的平均輪廓系數(shù)。結(jié)果顯示，當(dāng)分類數(shù)為5時(shí)，平均輪廓系數(shù)達(dá)到最大值0.57(圖6)，因此，將研究區(qū)內(nèi)數(shù)據(jù)分為5類較為合理，使用Turkey多重比較法檢驗(yàn)各類別間差異的顯著性，結(jié)果顯示各類別常量元素含量具有顯著性差異(P<0.05)，各類別樣品數(shù)及元素平均含量見(jiàn)表3。

圖6 不同分類數(shù)及對(duì)應(yīng)的平均輪廓系數(shù)Fig.6 Mean silhouette coefficients of differentclassification numbers

表3 各類別樣品常量元素平均值統(tǒng)計(jì)Table 3 Statistical table of mean values of major elements in various classifications %

為研究分類結(jié)果的實(shí)際內(nèi)涵及應(yīng)用價(jià)值，繪制了樣品類別的空間分布(圖7)。

圖7 樣品類別空間分布Fig.7 Spatial distribution map of sample categories

第Ⅰ類樣品數(shù)最少，僅為18個(gè)，但其元素含量特征明顯，SiO2含量最低，平均值僅56.15%，而Fe2O3、MgO、CaO含量最高(表3)，常量元素含量特征與基性巖較一致。第Ⅰ類樣品空間上位于五大連池玄武巖區(qū)，土壤類型為火山灰土。

第Ⅱ類樣品與第Ⅲ類樣品元素含量特征較為類似，其樣品數(shù)量也占絕對(duì)優(yōu)勢(shì)，分別占全區(qū)的28.5%和43.8%。第Ⅱ類樣品的SiO2含量(平均值64.66%)略高于第Ⅲ類樣品(平均值63.07%)，反映其相對(duì)含有更高的砂質(zhì)成分。在空間上，第Ⅱ類樣品主要分布在小興安嶺西麓丘陵地區(qū)，海拔普遍大于200 m(圖1)，該地區(qū)在第四紀(jì)時(shí)期以隆升為主，地形切割較劇烈，下部的砂層和砂礫石層多出露近地表，土壤類型以暗棕壤為主。

第Ⅲ類樣品數(shù)量最多，占全區(qū)的43.8%，其常量元素特征與典型黑土的成土母質(zhì)——黃土狀亞黏土[34]基本一致，廣泛分布在松遼平原東北部的高平原區(qū)，土壤類型多為黑土和草甸土。

第Ⅳ類樣品僅778件，常量元素表現(xiàn)為含有較豐富的CaO和MgO，其平均含量分別達(dá)5.04%和1.7%，代表土壤中碳酸鹽礦物含量較高?？臻g上沿個(gè)別河流分布，并均分布在隆起抬升的陡坡一側(cè)，與黑鈣土的分布吻合度較高。

第Ⅴ類樣品SiO2、K2O、Na2O含量較高，說(shuō)明土壤中含有較多的石英和長(zhǎng)石類礦物。在空間上，第Ⅴ類樣品主要分布在松花江南部各河流兩岸及一級(jí)階地。

3.4 地球化學(xué)分類結(jié)果與地質(zhì)單元的關(guān)系

土壤中元素地球化學(xué)特征受成土母質(zhì)、氣候、地形、生物、時(shí)間等因素綜合控制，而以往研究表明，土壤常量元素特征與成土母質(zhì)的關(guān)系最為密切[32,35]，土壤常量元素的空間分布與成土母質(zhì)和地質(zhì)單元往往具有較高的吻合度[36]。

本次聚類結(jié)果中，除第Ⅰ類樣品與玄武巖區(qū)完全對(duì)應(yīng)外，其他4類樣品在空間上與地質(zhì)單元的對(duì)應(yīng)性并不明顯(圖7、圖8)，其統(tǒng)計(jì)學(xué)結(jié)果也得出同樣結(jié)論(圖9)。原因可能包括兩個(gè)方面：一方面，地質(zhì)圖的比例尺與土壤地球化學(xué)調(diào)查尺度不匹配，在以往開(kāi)展的基礎(chǔ)地質(zhì)調(diào)查工作中，第四系廣泛覆蓋的平原區(qū)并不是調(diào)查重點(diǎn)，東北平原區(qū)目前編制的第四系地質(zhì)圖僅為1∶100萬(wàn)比例尺，精度較低；另一方面，同一地質(zhì)單元的地表巖性差異顯著，第四紀(jì)地質(zhì)單元通常由多個(gè)沉積地層構(gòu)成，如本區(qū)沖湖積物通常上部為亞黏土，下部則為砂層、砂礫層，受地表侵蝕、沉積、風(fēng)化等作用影響，同一地質(zhì)單元在不同地區(qū)出露巖性有較大差異，因此，地質(zhì)單元與地球化學(xué)的對(duì)應(yīng)關(guān)系并不完全一致[37-38]，也是造成二者對(duì)應(yīng)性低的一個(gè)原因。

圖8 典型黑土區(qū)第四紀(jì)地質(zhì)Fig.8 Quaternary geological map of typical black soil area

圖9 不同類別所包含地質(zhì)單元的累積頻率Fig.9 Cumulative frequency map of geologicalunits in different classifications

由于區(qū)域多目標(biāo)地球化學(xué)調(diào)查為1∶25萬(wàn)比例尺，采樣密度大，因此地球化學(xué)分類結(jié)果比地質(zhì)圖能更清晰地反映成土母質(zhì)的空間分布情況。綜合地理、地質(zhì)和土壤特征，可將5個(gè)類別分別命名為：Ⅰ—五大連池玄武巖母質(zhì)發(fā)育的火山灰土；Ⅱ—小興安嶺山前丘陵沖洪積物母質(zhì)發(fā)育的暗棕壤；Ⅲ—松遼平原東部高平原低鈣黃土狀土母質(zhì)發(fā)育的黑土；Ⅳ—松遼平原東部高平原邊緣富鈣沖湖積母質(zhì)上發(fā)育的黑鈣土；Ⅴ—松遼平原東南部沖積物母質(zhì)上發(fā)育的黑土。該結(jié)果較客觀地反映了典型黑土區(qū)的成土母質(zhì)分布情況。

3.5 地球化學(xué)分類結(jié)果反映的黑土生態(tài)環(huán)境問(wèn)題

地球化學(xué)分類結(jié)果表明，松花江南部土壤中砂質(zhì)含量較北部明顯偏高(圖7)，原因一方面是由于南部成土母質(zhì)以河流沖積物、沖洪積物和殘坡積物為主(圖8)，普遍含有較高的砂質(zhì)成分，但另一方面可能也反映出局部地區(qū)的黑土存在沙化現(xiàn)象。近期研究表明，中國(guó)東北部地區(qū)已成為沙塵暴的多發(fā)區(qū)，尤其是松花江、嫩江流域的局部草地退化和沙化地區(qū)[39]。受季風(fēng)和沙塵暴作用影響，每年都有大量細(xì)粒沙塵由西向東吹至黑土區(qū)。以哈爾濱地區(qū)為例，現(xiàn)代塵暴中干沉降和濕沉降粉塵的REE模式和Sr-Nd同位素組成都表明，沙塵來(lái)源于科爾沁沙地和渾善達(dá)克沙地[40-41]，這些沙塵的加入會(huì)增加黑土中砂質(zhì)含量，從而造成土壤沙化。在本次研究結(jié)果中，松花江南部沖湖積母質(zhì)上發(fā)育的黑土中目前也具有較高的SiO2含量，指示這些地區(qū)土壤可能已發(fā)生沙化，未來(lái)應(yīng)作為黑土區(qū)的生態(tài)問(wèn)題加以關(guān)注。

4 結(jié)論

1)典型黑土區(qū)常量元素SiO2、Al2O3和Fe2O3含量占80%以上，具有明顯的硅鋁土特點(diǎn)，表層土壤和深層土壤的常量元素比以0.9～1.1為主，說(shuō)明土壤常量元素主要受控于成土母質(zhì)。

2)利用主成分聚類法將黑土區(qū)土壤樣品劃分為5類最為合理，各類別地球化學(xué)特征差異顯著，其空間分布能較清晰地反映出成土母質(zhì)的空間分布情況。根據(jù)地理、地質(zhì)及土壤特征分別命名為：I—五大連池玄武巖母質(zhì)發(fā)育的火山灰土；Ⅱ—小興安嶺山前丘陵沖洪積物母質(zhì)發(fā)育的暗棕壤；Ⅲ—松遼平原東部高平原低鈣黃土狀土母質(zhì)發(fā)育的黑土；Ⅳ—松遼平原東部高平原邊緣富鈣沖湖積母質(zhì)上發(fā)育的黑鈣土；Ⅴ—松遼平原東南部沖積物母質(zhì)上發(fā)育的黑土。

3)松花江南部沖湖積母質(zhì)發(fā)育的黑土局部具有較高的SiO2含量，指示該地區(qū)黑土已出現(xiàn)沙化生態(tài)問(wèn)題，在黑土地保護(hù)中應(yīng)給予關(guān)注。