基于地統計學的土壤污染源解析模型差異對比

王彬，姜坤，師華定，徐嘉禮，吳海波，羅春暉

（1.浙江益壤環保科技有限公司，浙江 紹興 312000；2.紹興文理學院，浙江 紹興 312000；3.生態環境部土壤與農業農村生態環境監管技術中心，北京 100012；4.中國環境科學研究院土壤與固體廢物環境研究所，北京 100012）

土壤是地球系統的組成部分，更是人類及動植物賴以生存且不可替代的自然資源[1-2]。隨著工業化、農業化的不斷推進，我國土壤重金屬污染形勢逐漸加重，各個區域的污染程度不均，污染成因也不盡相同[3-8]，污染源解析顯得尤為重要。近年來，國內外學者在土壤重金屬源解析方面進行了大量研究[9-11]，主要采用源排放清單法、化學質量平衡模型、多元統計模型、先進數學算法以及空間分析法。其中，源排放清單法是基于污染源重金屬投入通量直接建立污染源清單數據庫，例如SHI等[12]匯總了我國農田土壤重金屬的排放清單；化學質量平衡模型是基于污染源和土壤重金屬的元素組成，根據質量守恒定律建立平衡模型，例如LIU等[13]、CHEN等[14]分別利用元素比值法、同位素比值法對研究區進行了污染源解析；多元統計模型是基于土壤重金屬元素組成的源解析方法，例如杜展鵬等[15]采用PCAAPCS-MLR模型得到了污染源對滇池草海和外海水質的貢獻程度，柴磊等[16]基于PMF模型對蘭州耕地土壤重金屬污染源進行了解析，劉玲玲等[17]采用UNMIX模型對北京城區公園土壤重金屬污染源進行了解析；先進數學算法是基于土壤重金屬元素組成和其他變量的解析方法，如HU等[18]利用條件推斷樹和有限混合模型識別了珠江三角地區土壤重金屬污染來源和作用機制，宋志廷等[19]利用隨機森林模型構建了天津武清地區表層土壤重金屬源匯的量化關系；空間分析法是基于土壤重金屬元素組成和空間關系的解析方法，如任加國等[20]利用空間關聯指數法對沱江上源支流土壤重金屬的空間變異及其相關性進行了解析。

以上關于源解析的研究方法較為單一，且局限性較大，如源排放清單法需要大量的數據清單，同位素比值法只能針對特定污染物（Cd、Hg、Pb）進行溯源，先進數學算法無法判斷解析結果對數據空間變異性的解釋程度等，且目前在流域尺度上基于地統計學方法將多種源解析模型進行對比分析的研究極少，而松陰溪流域是典型的灌溉水型河流，其周邊土地利用方式包括了耕地、園地以及成林，重工業也較為發達，適用于土壤污染源解析研究，因此本研究以該流域為研究對象，將多元統計模型和空間分析方法相結合，對研究區土壤重金屬污染類型及其貢獻率進行深入剖析，并將不同模型的解析結果進行對比分析，旨在為區域性土壤污染源解析工作提供方法借鑒和參考。

1 材料與方法

1.1 研究區概況

松陰溪是甌江上游主要支流，位于浙江省境西南部，源出遂昌縣垵口鄉北園附近，先自南而北，后折向東南，流經松陽縣，在麗水市大港頭入甌江中游大溪。松陰溪全長114 km，河寬約100 m，流域面積2 055 km2，灌溉著浙南的松古盆地。流域內有耕地約10 666.7 hm2，人口約22萬人，占松陽縣總人口的90%以上。該流域周邊道路縱橫交錯，分布著大量的企業，尤其是制造業和采礦業居多，且附近居民長期使用流域內的水資源進行農田灌溉，對園地進行大量的農藥噴灑，其周邊土壤環境污染問題日趨嚴重。

1.2 樣品的采集、制備與測定

利用ArcGIS10.6軟件對研究區表層土壤點位進行布設，所有點位均布設在農用地（耕地、林地、園地），原則上按照每2 500 m×2 500 m網格布設一個點位，以網格中心點作為采樣點，根據周邊具體情況適當調整。共布設62個表層土壤點位，對所有點位做2 km緩沖區后劃定影響范圍，并進行適當調整，如圖1所示。根據布設點位的經緯度，利用GPS對每個計劃采樣點進行精確定位，誤差不超過30 m，采樣人員到達計劃采樣點位后，需觀察其是否符合土壤采樣的代表性要求，在允許范圍內優選采樣點，位移距離不超過100 m，且低洼地、陡坡地、住宅、道路、溝渠附近等不布設點位。

圖1 研究區采樣點位、企業分布及土地利用情況Figure 1 Sampling points,enterprise distribution and land use in the study area

樣品采集以確定點位為中心劃定采樣區域，一般為50 m×50 m，采用蛇形法進行混合取樣，采樣前先用鐵鏟切割一個大于取土量的25 cm深的土方，然后用木鏟去除鐵鏟接觸面后裝入樣品袋，垂直挖土，保持采樣量上下一致，之后將采集完的土壤樣品置于陰涼處保存。樣品制備時，首先取適量土壤樣品平鋪在玻璃板上，避免陽光直射且室內溫度不超過40℃，自然風干，去除雜質，用木棒將土塊粉碎后過60目尼龍篩，再過100目尼龍篩，混合均勻后待測。采用原子熒光法和等離子體質譜法（ICP-MS）測定重金屬含量，分析方法精密度和準確度采用國家土壤標準物質GSS-15和室內平行樣品進行質量控制。

1.3 研究方法

1.3.1 PCA-APCS-MLR模型

PCA-APCS-MLR（主成分-絕對主成分-多元線性回歸）模型主要是通過降維充分解讀樣本數據所攜帶信息的差別，最后根據多元線性回歸系數來計算各類污染源的貢獻率[21-22]。首先需要計算樣本濃度的主成分得分，其公式為：

式中：(Az)k為主成分的得分值，i為樣本元素所在列號，wij為第j主成分的因子矩陣系數，zk為k觀測點污染物濃度標準化后的值，ck為k觀測點污染物濃度，c為污染物濃度的平均值，σ為污染物濃度的標準差，由于(Az)k為標準化的值，不能直接用于計算主成分的原始貢獻，必須將標準化的因子得分轉化為非標準化的絕對因子得分，才能用于主成分對污染指標的貢獻分析，其計算公式為：

式中：APCSjk為絕對主成分得分值，(Az)jk為第j主成分的得分值，(A0)j為0值下主成分得分值，i為樣本元素所在列號，Sij為第j主成分的因子矩陣系數，(Z0)i為觀測點污染物濃度0時標準化后的值，ci為污染物濃度的平均值，σi為污染物濃度的標準差。

以污染物濃度作為因變量，以絕對主成分得分為自變量，建立多元線性回歸方程，并利用回歸系數計算污染源的貢獻率，其計算公式為：

式中：m代表某種污染源類型，aim為回歸系數，bi為常數項，PCim為針對元素i、污染源m的貢獻率，aim×為元素i的所有樣本絕對主成分因子得分平均值。

該模型可通過多元線性回歸方程對污染物含量進行反預測，從而得到每個樣點對應的預測值。一般采用標準誤差和決定系數來分析模型精度，計算公式為：

式中，SE為標準誤差，R為決定系數，n為樣本數量，xi和yi分別為污染物含量的實測值和預測值。SE越小，R越大，代表模型精度越高，反之越低。

1.3.2 PMF模型

PMF（正定矩陣因子分解）模型法是對因子分析模型的一種優化，該模型已被成功運用于土壤污染源解析方面，它不僅克服了因子分析中負貢獻的情況，且模型運行不需要屬于源成分譜，通過樣本不確定度來預估模型誤差，故其解析結果更有意義[23-26]。PMF模型將原始數據分解成因子成分矩陣（F）和因子貢獻矩陣（G）以及一個殘差矩陣（E），之后根據不確定性計算各組分的誤差，并通過最小二乘法計算排放源的貢獻度。PMF模型計算公式為：

PMF模型為得到最優的因子成分矩陣（F）和因子貢獻矩陣（G），從而定義了一個目標函數Q，使其趨于自由度最小。其目標函數計算公式為：

式中：Xij為第i個樣品中第j個元素的濃度，p為污染物來源數量，gik為源k對第i個樣品的貢獻，fkj為源k中第j個元素的濃度，eij為殘差矩陣即PMF模型中未能解釋樣品濃度矩陣Xij的部分，uij為第i個樣品中第j個元素濃度的不確定性大小。不確定度可以根據元素的實測濃度與相應的方法檢出限（MDL）的關系確定，如果元素濃度的實測值沒有達到或者剛好達到方法檢出限時，不確定度計算公式為：

當元素遠遠超出了相應的檢出限時，不確定度的計算公式為：

式中：σ為標準差，c為元素濃度。

該模型運行結果對各污染物的標準誤差和決定系數均進行了統計，并繪制了污染物含量預測值與實測值的受體擬合圖，用以分析該模型的精度。

1.3.3 UNMIX模型

UNMIX模型是一種適用于解決混合問題的數學模型，受體點被認為是未知源組分的線性組合，這些來源對每個樣本都貢獻了一個未知的數量[27-29]。該模型假定來源的構成和貢獻都是正的。利用特定物種選擇的濃度數據估計每個樣本的源個數、源組成以及源對該樣本的貢獻，即在m個源的N個土壤樣品中，n個被分析的物種中的某一物種j的濃度可以用公式表示：

式中：Cij代表第i個樣品中第j個物種的濃度；Ujk代表第j個物種在源k中的質量分數，表示源組成；Dik代表源k在第i個樣品中的總量，視為源貢獻率；S代表各個源組成的標準偏差。

該模型運行結果同PMF模型類似，可直接觀察到污染物含量預測值與實際值的擬合情況，從而進一步判斷其精度。

2 結果與討論

2.1 土壤重金屬含量統計

研究區土壤重金屬含量統計情況見表1。對62個土壤樣品的統計結果進行分析可知，土壤pH的平均值和中位值分別為4.93和4.95，說明研究區土壤整體上呈弱酸性。Cd、Hg、As、Pb、Cr、Cu、Zn、Ni 8種重金屬含量的平均值和中位值均未超過《土壤環境質量 農用地土壤污染風險管控標準（試行）》（GB 15618—2018）風險篩選值，說明研究區土壤重金屬含量大部分處于較低水平。僅Cd、As兩種重金屬含量的最大值超過了風險篩選值但低于風險管制值，說明這兩種重金屬對研究區造成了一定污染，其他重金屬污染較輕。從偏度來看，數據分布均為右偏度，但經過Log變換后均符合正態分布；從峰度來看，數據均比正態分布高聳狹窄，說明其數據更集中在平均值附近。

表1 松陰溪流域土壤重金屬含量描述性統計Table 1 Descriptive statistics of soil heavy metals content in Songyin creek

2.2 不同源解析模型結果

2.2.1 PCA-APCS-MLR模型結果

利用SPSS26.0軟件對重金屬進行因子分析前，需要先對樣本數據進行檢驗，其中KMO檢驗結果為0.68，滿足最低0.5的要求；Bartlett檢驗觀測值為311.42，顯著性為0.00，拒絕零假設，說明該樣本數據總體適合因子分析。對研究區樣本數據進行主成分提取，且選擇特征值大于1的因子作為主因子，最終得到3個主因子。為了更加突出主因子的載荷能力，使結果更易于解讀，采用最大方差法對主因子進行旋轉，其結果見表2。從表2可以看出，3個主因子特征值變化幅度較小，累計方差貢獻率為79.94%，反映了樣本攜帶的絕大部分的變異，未旋轉前3個主因子可解釋方差分別為40.83%、24.93%、14.18%，旋轉后變為36.05%、29.42%、14.47%，模型得到優化。從表2可以看出，旋轉后的成分矩陣可以更好地識別因子特征，第一主因子的特征元素為Cd、Pb、As、Cu、Zn，第二主因子的特征元素為Cr、Ni，第三主因子的特征元素為Hg。通過PCA-APCS-MLR模型進一步計算各污染因子對各個元素所有樣點的貢獻率，負貢獻率改為0，其計算結果見圖2。可以看出，大部分貢獻率較高的重金屬為主因子對應的特征元素，而第二主成分對Cu的貢獻率也較高。

圖2 PCA-APCS-MLR模型對各重金屬污染源的貢獻率Figure 2 Contribution rate of heavy metal pollution sources by PCA-APCS-MLR

表2 PCA-APCS-MLR模型主因子解及成分矩陣Table 2 Main factors solution and component matrix of PCA-APCS-MLR

2.2.2 PMF模型結果

利用EPA PMF5.0軟件將樣本數據及其不確定度（根據方法檢出限求得）輸入模型，觀察8種重金屬數據的信噪比，將比值＜0.5設置為“bad”，0.5～1.0設置為“weak”，比值＞1.0設置為“strong”。利用旋轉工具觀察不同旋轉因子Fpeak下目標函數Q值的大小，多次迭代運行，發現旋轉因子Fpeak=-0.5時Q值最小，最終模型結果見圖3。從樣點來看，3種污染源的貢獻變化幅度均較大，說明污染源的區域性差異較為顯著；從元素來看，污染源1對Hg的貢獻較大，污染源2對Cd、As的貢獻較大，對Pb、Cu、Zn也有一定的貢獻，污染源3對Cr、Ni的貢獻較大。

圖3 PMF模型對各重金屬污染源的貢獻率Figure 3 Contribution rate of heavy metal pollution sources by PMF

2.2.3 UNMIX模型結果

利用EPA UNMIX6.0軟件對樣本數據進行模擬，首先對數據進行處理，利用Suggest Exclusion工具自動將濃度方差大于50%的數據進行篩選，結果顯示8種重金屬中僅Hg的濃度方差略高于50%，由于該元素為源解析重要因素，因此也可納入模型。UNMIX模型是利用幾何的“邊”對樣本數據進行分析，以各重金屬的濃度和作為基底，單個重金屬元素作為因變量，從而分析各個組分的邊緣診斷結果，通過運行該模型可以發現，除Cd、Hg外都具有明顯的上邊界或下邊界，通常認為大部分元素具有邊界時，模型基本是可靠的，且Cd、Hg兩種元素均無缺失值，各組分信噪比均符合要求，因此保留所有組分。

將8種重金屬元素的濃度均作為初始物種，將濃度和（total）作為模型輸入組分的總物種和標準物種，計算出符合要求的3源方案，結果顯示MinRsq=0.86，MinSig/Noise=2.83，均符合模型要求，說明該模型可以解釋86%的物種方差，源解析方案結果可靠。利用Excel進一步計算源貢獻率，如圖4所示。可以看出，污染源1對Cr、Ni的貢獻較大，其他元素相對較小；污染源2對Hg的貢獻較大，對Pb、Zn也有一定的貢獻；污染源3對Cd、As的貢獻較大，對Pb、Cu、Zn也有一定的貢獻。

圖4 UNMIX模型對各重金屬污染源的貢獻率Figure 4 Contribution rate of heavy metal pollution sources by UNMIX

2.3 污染源空間識別

為進一步對污染源進行識別，利用ArcGIS10.6對前文3種模型的樣點貢獻值進行地統計插值分析，如圖5所示。對于PCA-APCS-MLR模型來說，污染源1貢獻值較高的大部分區域分布的企業較少，道路情況一般，周邊無地表水，且無其他明顯特征，因此判別為“自然源”；污染源2貢獻值較高的區域主要分布著大量企業，主要包括采礦業和制造業，相關研究表明[30-32]，企業在冶煉加工、機械制造、廢氣排放時均會造成重金屬污染，因此判別為“工業源”；污染源3貢獻值較高的區域主要分布著復雜的道路（省道、縣道以及城鎮街道），且建有大型的停車場，由上文可知，該污染源的特征元素是Hg，根據以往研究表明[33-35]，交通因素對Hg的累積具有較大的作用，因此判別為“交通源”。對于PMF來說，污染源1貢獻值較高的區域與PCA-APCS-MLR模型的污染源3類似，且部分區域屬于農田和果園，當地村民利用附近未處理過的地表水進行灌溉，且相關研究表明[36-37]，果園采用大量的化學藥劑也會導致土壤重金屬的累積，因此判別為“交通-農藥-污灌源”；污染源2貢獻值較高的區域無明顯特征，因此判別為“自然源”；污染源3與PCAAPCS-MLR模型的污染源2較為類似，因此判別為“工業源”。對于UNMIX模型來說，污染源1貢獻值較高的區域分布著大量的企業，因此判定為“工業源”；污染源2貢獻值較高的區域與PMF模型的污染源1類似，因此判定為“交通-農藥-污灌源”；污染源3貢獻值較高的區域無明顯特征，因此判定為“自然源”。

圖5 各源解析模型中污染源貢獻率的空間分布Figure 5 Spatial distribution of pollution source contribution rate in each source analysis model

2.4 源解析模型結果對比

為了更好地對比源解析模型結果，首先需要對其模型精度進行對比。其中PCA-APCS-MLR模型通過多元線性回歸計算各重金屬元素的標準誤差和決定系數，而PMF模型和UNMIX模型在運行結果中可直接觀察預測值和實測值的擬合曲線，并給出了標準誤差和決定系數（相關系數的平方），統計結果見表3。從元素上來講，除Pb、Zn外，其他元素的標準誤差均較低，且大部分元素的決定系數在0.7以上，說明各模型整體擬合效果較好；從模型上來講，各模型的標準誤差和決定系數較為相似，其中PMF模型的標準誤差相對較小，可能是由于該模型在分析時添加了樣本的不確定性數據，從而降低了整體的預測誤差。

表3 各源解析模型對重金屬元素的擬合精度Table 3 Fitting accuracy of each source analytical model to heavy metal elements

結合地統計學方法對各模型的污染源類型及其貢獻進行分析，可以發現3種模型之間存在一定的差異性和相似性，從污染源類型來講，通過以上3種源解析模型都能識別出工業源與自然源，但利用PCAAPCS-MLR模型僅識別出交通源，而利用PMF模型、UNMIX模型卻識別出交通-農藥-污灌源，這可能是由于PCA-APCS-MLR模型缺乏敏感性，對一些貢獻度較小的污染源不能甄別處理，只對主要的源成分進行篩選和識別，而PMF模型、UNMIX模型對源成分比較敏感，且對于相似的源成分進行歸類，因此兩種模型沒有將交通源、農藥源以及污灌源進行區分。從污染源貢獻來講，PCA-APCS-MLR模型的源總體貢獻率差異較大，而PMF模型、UNMIX模型差異較小，可能是由于通過PCA-APCS-MLR模型分析未識別出農藥源和污灌源，且容易出現負貢獻，使其他源對As、Zn、Cr、Ni、Hg的貢獻度特別高，而PMF模型、UNMIX模型事先假設了受體點沒有負貢獻，且對于源識別較為全面，因此其源總體貢獻率的計算較為準確合理。總體來說，PCA-APCS-MLR模型在判別主要污染源類型時具有較大優勢，但容易出現負貢獻，貢獻率容易出現偏高或偏低的情況；而PMF模型、UNMIX模型難以區分成分相似的源類型，但在計算源貢獻率時具有一定的優勢，不會出現負值或偏高或偏低的情況。因此，在土壤污染源解析中，有時需要將不同的受體模型結合起來，更加精確地對源類型及其貢獻率進行綜合判別和計算，對研究區土壤污染源類型和總貢獻率進行重新整合和計算，其中源類型遵循細分原則，源貢獻以PMF模型、UNMIX模型為主，得到最終結果為自然源30.1%、工業源37.4%、交通源18.0%、農藥-污灌源14.5%。

3 結論

（1）研究區土壤Cd、Hg、As、Pb、Cr、Cu、Zn、Ni 8種重金屬含量的平均值和中位值均未超過《土壤環境質量 農用地土壤污染風險管控標準（試行）》（GB 15618—2018）風險篩選值，整體處于較低水平。

（2）通過以上3種受體模型結果可知，自然源是導致研究區土壤Cd、Pb、As、Cu、Zn污染的主要因素，總貢獻率為30.1%；Cr、Ni污染的主要來源是工業源，總貢獻率為37.4%；Hg污染主要與交通源、農藥-污灌源有關，總貢獻率分別為18.0%、14.5%。

（3）通過對比3種模型發現，PCA-APCS-MLR模型在判別主要的污染源類型時具有較大優勢，而PMF模型、UNMIX模型在計算源貢獻率時較為準確。