基于高斯混合模型的土壤環(huán)境質(zhì)量分區(qū)研究

劉佳斌，郜允兵，李永濤，徐卓揆，潘瑜春，楊晶，張婉秋，4

（1.華南農(nóng)業(yè)大學(xué)，廣州 510642；2.北京農(nóng)業(yè)信息技術(shù)研究中心，北京 100097；3.長沙理工大學(xué)交通運(yùn)輸工程學(xué)院，長沙 410114；4.中國礦業(yè)大學(xué)（北京），北京 100083）

土壤重金屬污染已經(jīng)成為世界性的環(huán)境問題[1]。2014 年原環(huán)境保護(hù)部和原國土資源部聯(lián)合發(fā)布的《全國土壤污染狀況調(diào)查公報(bào)》指出耕地土壤重金屬點(diǎn)位超標(biāo)率為19.4%，耕地土壤環(huán)境質(zhì)量堪憂。2016年5 月國務(wù)院發(fā)布了《土壤污染防治行動(dòng)計(jì)劃》（簡稱“土十條”），要求深入開展土壤重金屬污染調(diào)查，掌握土壤環(huán)境質(zhì)量狀況，實(shí)施農(nóng)用地分類管理，保障農(nóng)業(yè)生產(chǎn)環(huán)境安全[2]。我國一直關(guān)注土壤環(huán)境質(zhì)量采樣工作，“七五”以來，環(huán)境、農(nóng)業(yè)、國土等部門開展了多次土壤環(huán)境監(jiān)測調(diào)查工作，積累了大量的土壤環(huán)境質(zhì)量數(shù)據(jù)[3]?；谕寥啦蓸狱c(diǎn)開展環(huán)境質(zhì)量分區(qū)研究可為土壤污染的精準(zhǔn)識(shí)別、風(fēng)險(xiǎn)管控區(qū)的精準(zhǔn)劃定、管控措施精準(zhǔn)制定提供輔助決策依據(jù)，是當(dāng)前土壤環(huán)境質(zhì)量評(píng)價(jià)技術(shù)研究的熱點(diǎn)，具有重要的學(xué)術(shù)價(jià)值和現(xiàn)實(shí)意義[4-5]。

土壤環(huán)境質(zhì)量受自然和人為因素影響，影響因素（成土母質(zhì)、土壤類型、成土過程、土壤質(zhì)地、河流、道路、重點(diǎn)企業(yè)等）與土壤環(huán)境質(zhì)量空間分異并不一致，直接采用環(huán)境影響因素進(jìn)行中小尺度的環(huán)境質(zhì)量分區(qū)效果不佳[6-7]。2002 年，美國環(huán)境學(xué)家DORAN 等[8]提出了土壤環(huán)境質(zhì)量的概念，揭示了固有環(huán)境條件、環(huán)境質(zhì)量指標(biāo)與環(huán)境質(zhì)量空間格局形成的關(guān)系。KARLEN 等[9]從土壤使用功能入手分析了微觀、宏觀尺度上導(dǎo)致土壤環(huán)境質(zhì)量變化影響因素。郭書海等[10]以土壤環(huán)境質(zhì)量影響因素為基礎(chǔ)，采用影響因素空間疊加分析法進(jìn)行了國家尺度的土壤環(huán)境質(zhì)量綜合區(qū)劃研究。

土壤采樣是對(duì)區(qū)域環(huán)境質(zhì)量的抽樣，其樣點(diǎn)蘊(yùn)含且反映了區(qū)域環(huán)境質(zhì)量的情況，基于采樣點(diǎn)的環(huán)境質(zhì)量分區(qū)，是區(qū)域環(huán)境質(zhì)量分類后在空間域上的映射，可以通過對(duì)采樣點(diǎn)的屬性/空間聚類來實(shí)現(xiàn)。常用的點(diǎn)模式聚類方法主要有基于劃分、基于密度、基于模型、基于網(wǎng)格、基于層次等方法，在使用時(shí)也存在一定局限性[11-14]。局部Moran’s I 指數(shù)常用于基于屬性相關(guān)性的點(diǎn)模式空間聚類，可以對(duì)樣點(diǎn)觀測值的高高（高值點(diǎn)聚集）、高低（高值點(diǎn)被低值點(diǎn)包圍，表現(xiàn)為高值異常點(diǎn)）、低高（低值點(diǎn)被高值點(diǎn)包圍，表現(xiàn)為低值異常點(diǎn)）、低低（低值點(diǎn)聚集）等聚集模式和異常鄰接模式識(shí)別[15]，用于對(duì)污染指標(biāo)的空間聚集和異常區(qū)進(jìn)行刻畫，但其聚類或者異常的空間模式僅是基于樣點(diǎn)均值相對(duì)高低的反映，與人們對(duì)真實(shí)環(huán)境污染閾值絕對(duì)量的認(rèn)知仍有區(qū)別[16]。K-均值家族聚類以類與類之間的空間距離之和最小為約束條件，單靠空間距離不能真正體現(xiàn)人為、自然影響綜合作用環(huán)境質(zhì)量的分異性和相似性[17]。土壤環(huán)境質(zhì)量可以看作自然環(huán)境因素、人為活動(dòng)因素等多個(gè)驅(qū)動(dòng)力的單獨(dú)或者交互作用，將環(huán)境影響因素對(duì)土壤指標(biāo)的作用近似看作多個(gè)高斯過程或者其他隨機(jī)過程，則基于采樣點(diǎn)的環(huán)境質(zhì)量分類分區(qū)可以采用高斯混合模型進(jìn)行聚類分析。徐燁等[18]運(yùn)用高斯混合模型（Gaussian mixed model，GMM）聚類方法將采樣點(diǎn)總體分解為多個(gè)不同環(huán)境效應(yīng)作用的高斯簇集合，探索了土壤采樣點(diǎn)類型與企業(yè)空間分布間的關(guān)聯(lián)。但該聚類算法的目標(biāo)函數(shù)初始分類參數(shù)設(shè)定存在一定主觀性，容易造成聚類結(jié)果的過擬合或欠擬合[19-21]。

常用的土壤環(huán)境質(zhì)量評(píng)價(jià)方法有單項(xiàng)污染指數(shù)法、內(nèi)梅羅指數(shù)法、污染負(fù)荷指數(shù)法、環(huán)境風(fēng)險(xiǎn)指數(shù)法、地累積指數(shù)法、潛在生態(tài)危害指數(shù)法、灰色聚類法、健康風(fēng)險(xiǎn)評(píng)價(jià)方法等。上述這些評(píng)價(jià)方法被單獨(dú)或者聯(lián)合應(yīng)用于農(nóng)田土壤污染評(píng)價(jià)，前人利用單項(xiàng)污染評(píng)價(jià)指數(shù)、內(nèi)梅羅綜合污染指數(shù)、地積累指數(shù)、潛在生態(tài)危害指數(shù)等方法[22-24]開展評(píng)價(jià)工作，并評(píng)估了農(nóng)田土壤污染的程度，計(jì)算了各種風(fēng)險(xiǎn)指數(shù)，但這些方法反映的是土壤環(huán)境質(zhì)量某方面或者對(duì)污染指標(biāo)加權(quán)綜合，不一定適合環(huán)境質(zhì)量分區(qū)識(shí)別。分析研究區(qū)土壤環(huán)境質(zhì)量指標(biāo)的最小集，也是一種環(huán)境質(zhì)量的表征方法?；诖?，本文提出基于采樣點(diǎn)與輔助變量的區(qū)域土壤環(huán)境質(zhì)量分區(qū)方法，采用主要影響因素對(duì)土壤環(huán)境質(zhì)量預(yù)分類，通過主成分分析獲得的特征向量集合表征土壤環(huán)境質(zhì)量，以此構(gòu)建高斯混合模型對(duì)土壤采樣點(diǎn)聚類，避免采樣點(diǎn)分類初始化的隨機(jī)性導(dǎo)致聚類結(jié)果陷入局部最優(yōu)解，并以北京市順義區(qū)土壤采樣數(shù)據(jù)為例驗(yàn)證方法的有效性。

1 材料與方法

1.1 研究區(qū)概況

順義區(qū)位于北京市東北部，位于40°00′～40°18′N、116°28′～116°59′E。區(qū)內(nèi)耕地面積為404.05 km2，占土地總面積的39.56%，為北京市傳統(tǒng)農(nóng)業(yè)區(qū)。主要種植作物為小麥、玉米、蔬菜和水果。農(nóng)業(yè)用地以糧食生產(chǎn)為主，種植模式為冬小麥-夏玉米。土壤樣品采集于2009 年8 月至11 月，樣本點(diǎn)的布局和數(shù)量根據(jù)農(nóng)用地的利用方式和面積進(jìn)行確定，采樣方法為4 分法混合，包含4 個(gè)子樣，共采集耕層（0～20 cm）土壤樣品412 個(gè)，使用GPS 定位記錄樣點(diǎn)中心位置，采樣點(diǎn)主要分布于糧田、菜地、果園、設(shè)施農(nóng)業(yè)用地等，其中后沙峪、南法信、天竺、馬坡、仁和等區(qū)域未進(jìn)行采樣，為非調(diào)查區(qū)。根據(jù)采樣點(diǎn)的土地利用方式，分為菜地、果園、荒地、林地、苗圃、設(shè)施農(nóng)業(yè)用地和水澆地，采樣點(diǎn)分布如圖1 所示。所有樣品在室內(nèi)自然風(fēng)干，碾壓磨碎后過100 目尼龍網(wǎng)篩，土壤樣品按照EPA Methods 3050B 的方法使用HNO3、HCl 和H2O2混合溶液消解。消解液中的As、Cd、Cu、Pb和Zn使用等離子發(fā)射光譜（ICP-OES，Thermo iCAP 6300，USA）測定，使用HNO3∶HCl（1∶1）在100 ℃消解，消解液使用原子熒光光譜（AFS，Titan AFS 830，China）測定Hg 的含量。每個(gè)樣品重復(fù)測3 次，樣品處理和分析的質(zhì)量控制采用分析土壤標(biāo)準(zhǔn)參考物質(zhì)GSS-1 和GSS-4（國家標(biāo)準(zhǔn)物質(zhì)中心）。

成土母質(zhì)[25]、土壤類型、土壤質(zhì)地、土地利用方式、工礦企業(yè)排放、河流等是可能影響該區(qū)域土壤環(huán)境質(zhì)量空間差異的因素[26-27]。研究區(qū)土壤成土母質(zhì)均為第四紀(jì)黃土，不存在巖石巖性上的差異。區(qū)域內(nèi)土壤類型包含潮土、褐土、水稻土3類，土壤質(zhì)地分為輕壤質(zhì)、中壤質(zhì)、砂壤質(zhì)，其空間分布如圖2（a）、圖2（b）所示，土壤類型的空間連片性強(qiáng)于土壤質(zhì)地；區(qū)內(nèi)主要河流有懷河、潮白河、金雞河等，工礦企業(yè)主要分布在河流下游區(qū)域附近，其空間分布如圖2（c）、圖2（d）所示，河流、工礦企業(yè)按照輻射范圍劃定影響區(qū)，可以看出工礦企業(yè)聚集度具有一定的空間異質(zhì)性。

1.2 研究方法

基于高斯混合模型的土壤環(huán)境質(zhì)量分區(qū)研究主要包括3 個(gè)部分：首先，對(duì)土壤采樣點(diǎn)的重金屬指標(biāo)歸一化處理后進(jìn)行主成分分析，選取特征值大于1 的主成分共同反映區(qū)域土壤環(huán)境質(zhì)量綜合特征，并運(yùn)用地理探測器篩選影響土壤環(huán)境質(zhì)量綜合特征的主要因素，按主要影響因素類別對(duì)采樣點(diǎn)進(jìn)行預(yù)分類；其次，基于采樣點(diǎn)的預(yù)分類方案構(gòu)建高斯混合模型（Gaussian mixed model，GMM），以期望最大化算法（Expectation-maximum algorithm，EM 算法）對(duì)高斯混合模型中各子高斯函數(shù)的權(quán)重（α）、均值（μ）、標(biāo)準(zhǔn)差（σ）參數(shù)進(jìn)行迭代求解；最后，以采樣點(diǎn)構(gòu)建泰森多邊形，并將采樣點(diǎn)分類結(jié)果賦值到多邊形后進(jìn)行同類多邊形合并，形成環(huán)境質(zhì)量初步分區(qū)，結(jié)合自然與人為影響因素的范圍對(duì)分區(qū)結(jié)果進(jìn)行邊界調(diào)整。

1.2.1 主要影響因素分析

區(qū)域土壤重金屬受自然或者人為因素影響，重金屬在產(chǎn)生過程中往往存在一定伴生性[28]，在數(shù)值上表現(xiàn)為強(qiáng)弱不一的相關(guān)性，直接采用原始觀測值進(jìn)行聚類分析容易受到噪聲干擾。本文采用主成分分析法（PCA）對(duì)歸一化后的412個(gè)土壤樣品的重金屬指標(biāo)（Cd、As、Cr、Hg、Cu、Pb、Zn等）進(jìn)行特征降維，篩選前N個(gè)數(shù)值＞1的特征值（N≥1）及對(duì)應(yīng)特征向量，并對(duì)原始數(shù)據(jù)進(jìn)行特征向量投影變換，變換后得到特征向量組PC（PC1，PC2，…PCN），N個(gè)特征向量可表征區(qū)域環(huán)境質(zhì)量綜合特征，反映不同重金屬指標(biāo)的總體相似性和差異性。

采樣點(diǎn)的初始分類數(shù)對(duì)最終環(huán)境質(zhì)量分類有影響，容易導(dǎo)致分類結(jié)果陷入局部最優(yōu)解。在進(jìn)行高斯混合建模前有必要對(duì)區(qū)域環(huán)境質(zhì)量主要環(huán)境影響因素進(jìn)行識(shí)別，通過主要影響因素分類對(duì)采樣點(diǎn)進(jìn)行預(yù)分類。地理探測器[29]可以用于探測環(huán)境質(zhì)量特征變量PCN的空間分異性（也稱空間分層或者分類），即探測環(huán)境質(zhì)量影響因素Y（是土壤類型、土壤質(zhì)地、土地利用方式、工礦企業(yè)排放、河流等中的一個(gè)類別型變量）多大程度上解釋了環(huán)境質(zhì)量特征變量PCN的空間分異，統(tǒng)計(jì)量q通過F檢驗(yàn)，則說明該影響因素有較強(qiáng)解釋力。具體計(jì)算公式及步驟如下：

式中：h=1，2，…，L為環(huán)境質(zhì)量影響因素Yi（z1i，z2i，…，zLi）分類或環(huán)境質(zhì)量特征PC 的分層；Nh和N分別為層h和全區(qū)的采樣點(diǎn)數(shù)分別為相應(yīng)層h和全區(qū)的土壤環(huán)境質(zhì)量特征的方差。q?[0,1]，其q值越大表明影響因素Y對(duì)特征變量PCN的解釋力越強(qiáng)，此時(shí)環(huán)境影響因素Y的分類數(shù)可以作為高斯混合模型建模的預(yù)分類數(shù)K。q的顯著性可采用F非中心分布進(jìn)行檢驗(yàn)，見公式（2）。

式中：N為采樣點(diǎn)數(shù)；λ為非中心度。

對(duì)研究區(qū)域中各環(huán)境影響因素分層下的異質(zhì)性進(jìn)行比較，通過顯著性檢驗(yàn)的影響因素為采樣點(diǎn)的主要影響因素。由于每個(gè)采樣點(diǎn)空間位置上的土壤環(huán)境輔助變量因子屬性具有唯一性，因此通過空間對(duì)應(yīng)關(guān)系可將具有相同影響因素類型的采樣點(diǎn)分為同一類，作為土壤采樣點(diǎn)的預(yù)分類。該預(yù)分類方法基本可以保證各類采樣點(diǎn)之間呈現(xiàn)出分異規(guī)律，是一種合理的初始化方式。

1.2.2 采樣點(diǎn)的聚類方法

高斯混合模型是多個(gè)高斯分布函數(shù)的線性組合，理論上GMM可以擬合出任意類型的分布，通常用于解決同一集合下的數(shù)據(jù)包含多個(gè)不同的分布、同一類分布但參數(shù)不同、不同類型的分布（比如正態(tài)分布和伯努利分布）等情況。區(qū)域環(huán)境質(zhì)量受自然和人為因素綜合作用，其概率分布往往呈右偏態(tài)分布（圖3），若將自然或者人為因素影響的過程看作多個(gè)高斯過程或者其他隨機(jī)過程的綜合效應(yīng)，則可用高斯混合模型表示。高斯混合模型一般采用EM算法進(jìn)行參數(shù)的求解。本文以主要影響因素的K個(gè)分層，即主要影響因素的K個(gè)類型對(duì)土壤環(huán)境質(zhì)量綜合特征預(yù)分類（見1.2.1），在此基礎(chǔ)上，通過高斯混合模型對(duì)特征變量PC聚類建模，并使用期望最大化算法（EM算法）對(duì)模型中的參數(shù)進(jìn)行估計(jì)，具體步驟如下：

（1）高斯混合模型初始化

將K個(gè)高斯模型線性組合的高斯混合模型P作為區(qū)域環(huán)境質(zhì)量模擬，其表現(xiàn)形式為：

式中：a=(a1,a2,…,an),ai為第i個(gè)土壤采樣點(diǎn)上的環(huán)境綜合特征向量；K為子高斯模型個(gè)數(shù)；θ=（α1,α2,…,αK;θ1,θ2,…,θK）；αk為系數(shù)為第k個(gè)分模型高斯分布密度分別為第k個(gè)分布的期望和方差。

初始化時(shí)以模型變量γik表示分類結(jié)果，結(jié)合采樣點(diǎn)與第k類環(huán)境影響因素因子的對(duì)應(yīng)關(guān)系，預(yù)分類過程如下：

所得到的γik初始值由采樣點(diǎn)在土壤環(huán)境主導(dǎo)因素影響下的分類情況確定。則完整數(shù)據(jù)集合A=完整數(shù)據(jù)的對(duì)數(shù)似然函數(shù)為：

（2）參數(shù)估算

E步：計(jì)算分模型k對(duì)觀測數(shù)據(jù)ai的響應(yīng)度，計(jì)算方法如下：

其中θ是上一步迭代運(yùn)算求得的參數(shù)值。

M步：求Q對(duì)θ的極大值，即求新一輪迭代的模型參數(shù)，計(jì)算方法如下：

1.2.3 初始分類結(jié)果的分區(qū)

在實(shí)際土壤環(huán)境質(zhì)量管理中需要基于同類相近原則對(duì)土壤采樣點(diǎn)劃類歸區(qū)，形成具有等級(jí)關(guān)系的子區(qū)域[30]。本文在土壤采樣點(diǎn)聚類分析結(jié)果上，根據(jù)采樣點(diǎn)集的幾何位置生成相應(yīng)泰森多邊形（Voronoi圖），將具有共享邊（Rook 連接）且有相同聚類等級(jí)的泰森多邊形要素合并為初步分層要素，并賦予相應(yīng)的類型，實(shí)現(xiàn)土壤采樣點(diǎn)“點(diǎn)到面”的變換（泰森多邊形法）。針對(duì)破碎化現(xiàn)象嚴(yán)重的區(qū)域，依據(jù)相近相似性原則，進(jìn)一步結(jié)合河流、道路等邊界對(duì)分區(qū)范圍進(jìn)行邊界調(diào)整。

1.3 數(shù)據(jù)處理過程

本文數(shù)據(jù)處理方法及步驟如下：（1）運(yùn)用SPSS 25對(duì)Cd、As、Cr、Hg、Cu、Pb、Zn 7種土壤重金屬含量進(jìn)行統(tǒng)計(jì)分析、正態(tài)分布檢驗(yàn)和皮爾森（Pearson）相關(guān)性分析，探索數(shù)據(jù)的統(tǒng)計(jì)規(guī)律和各重金屬之間的相關(guān)性。（2）通過SPSS 對(duì)多種土壤重金屬采樣數(shù)據(jù)進(jìn)行主成分分析，提取表征土壤環(huán)境質(zhì)量的主成分綜合變量。（3）采用地理探測器軟件（www.geodetector.cn）分析土壤類型、土壤質(zhì)地、土地利用類型等環(huán)境影響輔助因素與主成分向量之間的相關(guān)關(guān)系。（4）基于MATLAB 2017a 平臺(tái)編寫基于先驗(yàn)環(huán)境影響主導(dǎo)因素預(yù)分類的高斯混合模型（GMM）算法，對(duì)代表土壤環(huán)境質(zhì)量綜合特征的主成分向量進(jìn)行聚類分析。（5）根據(jù)采樣點(diǎn)空間位置信息運(yùn)用ArcGIS 10.5 生成相應(yīng)范圍的泰森多邊形，并將互為鄰域且具有相同分類屬性的多邊形要素合并后，作為初步分層（分類結(jié)果）。（6）在ArcGIS 平臺(tái)中結(jié)合研究區(qū)河流、道路分布，對(duì)初步分區(qū)結(jié)果的邊界進(jìn)行調(diào)整。

2 結(jié)果與分析

2.1 采樣點(diǎn)探索性分析

與北京地區(qū)土壤重金屬背景值相比較[31]，除Pb外其他重金屬（類金屬）指標(biāo)具有不同程度的積累效應(yīng)。Cd、As、Hg、Cu、Pb 等重金屬指標(biāo)的變異系數(shù)均大于25%，達(dá)到了中等變異水平。除Cr 外其他土壤重金屬頻率分布圖呈現(xiàn)不同程度右偏態(tài)，表明研究區(qū)內(nèi)局部地區(qū)受工礦企業(yè)、河流灌溉等因素影響，Hg、Cd、As指標(biāo)的最大值已接近篩選值，表明在某些樣點(diǎn)上存在一定程度的風(fēng)險(xiǎn)。北京市土壤pH 均值為7.64[32]，總體上看，各重金屬含量均值、中位值低于《土壤環(huán)境質(zhì)量農(nóng)用地土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)（試行）》（GB 15618—2018）中pH 值大于7.5 時(shí)的篩選值，從土壤污染評(píng)價(jià)的角度看該區(qū)土壤總體質(zhì)量屬于優(yōu)先保護(hù)類。順義區(qū)土壤重金屬描述性統(tǒng)計(jì)如表1所示。

表1 區(qū)域土壤重金屬統(tǒng)計(jì)特征（mg·kg-1）Table1 Statistical characteristics of heavy metal content in regional soil（mg·kg-1）

2.2 采樣點(diǎn)預(yù)分類結(jié)果

基于Pearson 系數(shù)法分析各元素之間的相關(guān)性，結(jié)果如表2所示，其中Zn與Cu之間的相關(guān)性最強(qiáng)，相關(guān)系數(shù)高達(dá)0.476，說明該兩個(gè)指標(biāo)具有相對(duì)較高的信息重疊，As 與Cd、As 與Cu、Cd 與Cu、Cd 與Pb 之間具有較為顯著的相關(guān)性（r≥0.2，P＜0.05），說明其相互之間也存在部分信息重疊，存在一定程度的同位或伴生性，有相關(guān)性的土壤重金屬指標(biāo)可能受相同自然或人類活動(dòng)因素影響。

表2 土壤重金屬含量相關(guān)系數(shù)矩陣Table2 The correlation coefficient matrix of soil heavy metal content

本研究采用PCA 分析法對(duì)土壤采樣點(diǎn)重金屬指標(biāo)進(jìn)行降維，其結(jié)果如表3 所示。選取特征值大于1.0（特征值表示原始指標(biāo)投影到該特征向量后的方差大小，方差越大則保留的原始信息量也越大）的對(duì)應(yīng)特征向量PC1（特征值1.885）、PC2（特征值1.112），二者累積貢獻(xiàn)率達(dá)到了42.8%，從表3 中主成分的因子載荷系數(shù)可以看出主成分PC1 主要表征土壤中As、Pb、Zn的含量，主成分PC2主要反映土壤中Cd、Cr的含量，可用PC1、PC2 特征向量表征區(qū)域土壤環(huán)境綜合質(zhì)量特征。

表3 土壤重金屬的變量解釋與主成分矩陣Table3 Variable interpretation and principal component matrix of heavy metals in soil

本文在采用主成分降維后的特征向量PC1、PC2進(jìn)行采樣點(diǎn)預(yù)分組時(shí)，通過地理探測器分析土壤環(huán)境影響因素Y{土壤質(zhì)地、土壤類型、土地利用類型、河流影響、工礦企業(yè)影響范圍}對(duì)特征向量PC1、PC2 的解釋力，即分析土壤環(huán)境影響因素Y對(duì)特征向量PC1、PC2的q值的顯著性。地理探測器檢驗(yàn)各分組下樣點(diǎn)主成分間分異結(jié)果如表4 所示。主成分PC1 在土壤質(zhì)地（輕質(zhì)壤、中質(zhì)壤、沙質(zhì)壤，根據(jù)砂粒度與黏粒度分類）、河流影響區(qū)的分類下q統(tǒng)計(jì)量滿足F檢驗(yàn)的顯著性差異。其中，河流影響區(qū)以相關(guān)研究[33]為基礎(chǔ)，通過半徑1 km 的河流緩沖區(qū)獲得。主成分PC2在土壤質(zhì)地、土壤類型與土地利用類型的分類下q統(tǒng)計(jì)量具有顯著性差異。研究區(qū)內(nèi)土壤母質(zhì)類型差異較小，而土壤質(zhì)地的黏粒度越重則重金屬淋溶程度越低，同時(shí)土壤膠體的吸附能力也越強(qiáng)[34]，在一定程度上影響了土壤重金屬的遷移和擴(kuò)散能力。土壤環(huán)境質(zhì)量的特征向量PC1、PC2 在土壤質(zhì)地各類型中皆具有顯著性差異，表明土壤質(zhì)地對(duì)土壤環(huán)境質(zhì)量空間分異有較好解釋力。因此，采用土壤質(zhì)地的空間分層作為土壤采樣點(diǎn)預(yù)分類方案。

表4 不同影響因素分組下樣點(diǎn)主成分間分異檢驗(yàn)結(jié)果Table4 Test results of soil samples in different groups

2.3 土壤環(huán)境質(zhì)量分區(qū)

采用土壤質(zhì)地類型分區(qū)（輕壤質(zhì)、中壤質(zhì)、砂壤質(zhì)）將土壤采樣點(diǎn)初步分成3 類，作為高斯混合模型構(gòu)建中研究區(qū)土壤環(huán)境質(zhì)量主成分特征向量的初始類型，利用EM 算法對(duì)高斯混合模型中未知參數(shù)進(jìn)行迭代求解，直到各子數(shù)據(jù)集符合正態(tài)分布概率最大。EM 算法對(duì)高斯混合模型（也稱為GEM 聚類）進(jìn)行求解的終止條件為似然函數(shù)收斂。經(jīng)過32 次EM 迭代運(yùn)算后，土壤采樣點(diǎn)分成A、B、C 3 個(gè)種類，各類樣點(diǎn)數(shù)目分別為279、99、34。基于此，運(yùn)用泰森多邊形劃分法構(gòu)建包含土壤采樣點(diǎn)的土壤環(huán)境質(zhì)量初步分區(qū)，結(jié)果如圖5（a）所示。參考河流、道路對(duì)初步分區(qū)中細(xì)碎圖斑進(jìn)行適度調(diào)整，調(diào)整后的局部單元內(nèi)部有多個(gè)類型空間分層時(shí)，將其歸為面積占優(yōu)的空間分區(qū)類型，最終形成土壤環(huán)境質(zhì)量分區(qū)[圖5（b）]。

為了進(jìn)一步說明研究區(qū)環(huán)境質(zhì)量分區(qū)合理性，本文查閱了影響土壤環(huán)境質(zhì)量的文獻(xiàn)資料[35-36]（禽養(yǎng)殖業(yè)分布信息、河流道路空間分布數(shù)據(jù)、土壤類型數(shù)據(jù)、工礦企業(yè)空間分布等），并以呈現(xiàn)右偏態(tài)分布的土壤重金屬指標(biāo)Pb 在各分區(qū)中統(tǒng)計(jì)分布為例展開說明，對(duì)三類區(qū)的環(huán)境質(zhì)量特點(diǎn)及可能原因進(jìn)行分析。（1）A 類環(huán)境質(zhì)量區(qū)：區(qū)內(nèi)土壤重金屬指標(biāo)Pb 均值為14.06 mg·kg-1，方差為2.19 mg·kg-1，能通過正態(tài)分布檢驗(yàn)，均值比較平穩(wěn)。該區(qū)內(nèi)土壤類型以潮土為主，面積占整個(gè)A 區(qū)面積的74.27%，土壤質(zhì)地以輕壤質(zhì)為主，占該子區(qū)域面積的50.80%，且與A 類分區(qū)具有較強(qiáng)的一致性，根據(jù)相關(guān)研究，輕質(zhì)潮土上的腐殖質(zhì)積累過程較弱[37]，河流影響區(qū)占該子區(qū)域的24%，區(qū)內(nèi)人為因素影響較小，土壤環(huán)境質(zhì)量相對(duì)穩(wěn)定；（2）B類環(huán)境質(zhì)量區(qū)：區(qū)內(nèi)土壤重金屬指標(biāo)Pb 均值為21.70 mg·kg-1，方差為3.72 mg·kg-1，比較平穩(wěn)，能通過正態(tài)分布檢驗(yàn)，但均值略高于A 類區(qū)，含量波動(dòng)情況也略大于A類區(qū)。該區(qū)內(nèi)同樣以輕壤質(zhì)潮土為主，潮土面積占整個(gè)B 區(qū)面積的49.42%，輕壤質(zhì)占該子區(qū)域面積的55.40%，相比A 區(qū)略有減少，區(qū)內(nèi)水資源相對(duì)豐富，河流總數(shù)和水量較大，按河流灌溉可達(dá)性分析占整個(gè)區(qū)內(nèi)面積的27%，該區(qū)域以農(nóng)業(yè)種植為主，主要種植水果蔬菜，用地呈現(xiàn)出高生產(chǎn)資料投入與輪作頻率大的特征，同時(shí)區(qū)內(nèi)農(nóng)村居民地分布較多[38]，土壤中容易發(fā)生一定的重金屬累積；（3）C 類環(huán)境質(zhì)量區(qū)：C 類區(qū)的土壤重金屬Pb均值和方差分別為27.97 mg·kg-1與18.91 mg·kg-1，變化波動(dòng)較大，但也滿足正態(tài)分布檢驗(yàn)。該區(qū)也以輕壤質(zhì)潮土為主，潮土面積占整個(gè)C 區(qū)面積的65.17%，輕壤質(zhì)占該子區(qū)域面積的50.80%，該區(qū)內(nèi)工礦企業(yè)影響面積占該區(qū)49.4%，工礦企業(yè)用地面積相對(duì)較大，且子區(qū)域處于潮白河、箭桿河主要干道上，受河流影響面積占比高達(dá)48.1%，同時(shí)在2009 年前后禽類養(yǎng)殖業(yè)也較發(fā)達(dá)，糞便和飼料粉塵通過污水排放和大氣沉降對(duì)周邊產(chǎn)生影響，易造成重金屬聚集[39]。

3 討論

為了進(jìn)一步說明聚類方法對(duì)土壤環(huán)境質(zhì)量分區(qū)結(jié)果的影響，本文采用K-means 聚類分析、自組織特征映射神經(jīng)網(wǎng)絡(luò)（Self-organizing feature maps，SOFM）聚類進(jìn)行土壤樣點(diǎn)聚類方法的對(duì)比分析。具體聚類實(shí)驗(yàn)時(shí)，也采用上文提到的基于土壤質(zhì)地預(yù)分類的土壤樣點(diǎn)的PC1、PC2 特征向量作為聚類模型計(jì)算的輸入，分類數(shù)K統(tǒng)一設(shè)置為3。在K-means聚類時(shí)，選擇歐式距離為度量方式，直到?jīng)]有對(duì)象被重新分配給不同的聚類，整體誤差平方和局部最小為止；在SOFM聚類中，競爭層神經(jīng)元個(gè)數(shù)依次設(shè)置為3，采用六角形拓?fù)浣Y(jié)構(gòu)，訓(xùn)練次數(shù)為500，直到競爭型神經(jīng)網(wǎng)絡(luò)取得最優(yōu)解，聚類中心沒有再發(fā)生變化。初步聚類結(jié)果按A、B、C 3類標(biāo)記到土壤樣點(diǎn)，采用泰森多邊形法構(gòu)建土壤樣點(diǎn)初步分區(qū)，并依據(jù)河流、道路對(duì)細(xì)碎圖斑進(jìn)行分區(qū)邊界微調(diào)。K-means 聚類下的A、B、C 類型點(diǎn)位數(shù)分別為244、56、112，在調(diào)整過程中有23.27 km2的B 類、12.55 km2的C 類調(diào)入A 類型，調(diào)整面積占總面積的2.8%；SOFM聚類下的A、B、C類型點(diǎn)位數(shù)分別為151、134、127，在調(diào)整過程中有16.13 km2的B類、46.88 km2的C 類調(diào)入A 類型，調(diào)整面積占總面積的7.73%，分區(qū)結(jié)果如圖6（a）、圖6（b）所示。

為比較3 種聚類分區(qū)結(jié)果的差異性，分別對(duì)各類型分區(qū)中的主成分進(jìn)行統(tǒng)計(jì)分析，不同分區(qū)下主成分特征向量統(tǒng)計(jì)特征如表5 所示，通過比較各成分特征向量在不同子分區(qū)中的均值差異，其中以GEM 聚類分區(qū)中的差異最為明顯；且在本文的GEM 聚類分區(qū)中，成分特征向量在A、B、C 3種分區(qū)下均能滿足顯著性檢驗(yàn)，能揭示受相同環(huán)境影響因素作用的可能性，而在SOFM 聚類分區(qū)中，主成分PC1 在A、C 類區(qū)中都不能通過顯著性檢驗(yàn)，在K-means 聚類分區(qū)中，主成分PC2在B類區(qū)中不能通過顯著性檢驗(yàn)。

表5 不同分區(qū)下主成分特征向量統(tǒng)計(jì)特征Table5 Statistical characteristics of principal component eigenvectors in different partitions

為進(jìn)一步比較3 種聚類方法的分區(qū)結(jié)果的差異性，本文采用地理探測器探測分區(qū)結(jié)果分析土壤環(huán)境綜合質(zhì)量的PC1、PC2 的解釋力。采用GMM 分區(qū)結(jié)果、SOFM分區(qū)、K-means分區(qū)中的類型數(shù)作為地理探測器中自變量PC 的分層數(shù)，PC1、PC2作為自變量，各分區(qū)下樣點(diǎn)綜合質(zhì)量分異性結(jié)果如表6 所示，其中主成分PC1 在GEM 聚類分區(qū)與SOFM 聚類分區(qū)下具有顯著性差異，主成分PC2 僅在GEM 聚類分區(qū)方法下具有顯著性差異，而K-means聚類方法下的主成分均無顯著性差異。其可能原因是：SOFM 聚類是基于神經(jīng)網(wǎng)絡(luò)訓(xùn)練的分類，在各權(quán)向量分別向各個(gè)聚類模式群的中心位置靠攏時(shí)，容易受神經(jīng)元排列的拓?fù)浣Y(jié)構(gòu)影響，導(dǎo)致權(quán)重分配不均[40]；針對(duì)土壤采樣點(diǎn)含量通常呈現(xiàn)右偏態(tài)分布的特性，K-means 聚類分析是以類與類之間的距離之和最小為約束的分類方法，只對(duì)各類的緊湊性進(jìn)行約束，卻不是針對(duì)重金屬含量的分布特征進(jìn)行劃分[41]；本研究中所采用的GEM 聚類方法是以多個(gè)高斯函數(shù)對(duì)右偏態(tài)分布的土壤環(huán)境質(zhì)量進(jìn)行描述，假設(shè)每個(gè)正態(tài)分布函數(shù)代表不同分類的重金屬含量，相同分類的采樣點(diǎn)具有同一數(shù)學(xué)期望和方差，將采樣點(diǎn)分成受不同環(huán)境因素作用的多個(gè)正態(tài)分布類型，GEM 聚類分析讓異類之間的差異性更大，同類之間的相似性更大，因而相比其他分層分區(qū)方法有著更好的空間分層異質(zhì)性，使所得的分區(qū)結(jié)果與土壤環(huán)境質(zhì)量空間分布具有高度一致性。

表6 不同分區(qū)方法下主成分間分異檢驗(yàn)結(jié)果Table6 Test results of differentiation between principal components under different zoning methods

4 結(jié)論

（1）土壤重金屬含量是環(huán)境影響因素綜合作用的結(jié)果，基于高斯混合模型的聚類方法可用于區(qū)域環(huán)境質(zhì)量分區(qū)?；谥饕绊懸蛩氐耐寥啦蓸狱c(diǎn)預(yù)分類策略，避免了土壤采樣點(diǎn)初始化分類的盲目性，有助于聚類結(jié)果的快速收斂，避免聚類結(jié)果陷入局部最優(yōu)，環(huán)境質(zhì)量分區(qū)的結(jié)果更接近客觀認(rèn)知。該土壤環(huán)境質(zhì)量分區(qū)方法兼顧了采樣點(diǎn)屬性特征與環(huán)境輔助變量信息，分區(qū)方法在應(yīng)用中具有可推廣性，可為全國監(jiān)測點(diǎn)的環(huán)境質(zhì)量分區(qū)提供參考。

（2）以聚類分析的思想解析土壤環(huán)境質(zhì)量空間分布具有良好的效果，基于高斯混合模型的聚類分區(qū)充分考慮工礦企業(yè)、河流、道路、土壤質(zhì)地、土壤類型等不同影響的作用范圍和尺度特征，可準(zhǔn)確客觀反映土壤環(huán)境質(zhì)量分布的情況，是了解污染狀況、過程、格局的手段和方法。不同類型空間分類分級(jí)是進(jìn)行精準(zhǔn)管理的前提，準(zhǔn)確的先驗(yàn)分區(qū)可用于進(jìn)一步指導(dǎo)布點(diǎn)監(jiān)測調(diào)查。