基于貝葉斯最大熵和輔助信息的土壤重金屬含量空間預(yù)測(cè)

費(fèi)徐峰，任周橋，樓昭涵，肖銳，呂曉男*

（1.浙江省農(nóng)業(yè)科學(xué)院數(shù)字農(nóng)業(yè)研究所，杭州 310021；2.農(nóng)業(yè)農(nóng)村部農(nóng)產(chǎn)品信息溯源重點(diǎn)實(shí)驗(yàn)室，杭州 310021；3.浙江大學(xué)海洋學(xué)院，浙江 舟山 316021；4.武漢大學(xué)遙感與信息工程學(xué)院，武漢 430072）

土壤是人類生存和發(fā)展的物質(zhì)基礎(chǔ)，是糧食安全和人類健康的保證。近幾十年來(lái)，隨著快速的城市化和工業(yè)化進(jìn)程，包括重金屬在內(nèi)的土壤污染物增加已成為社會(huì)發(fā)展和人類健康的重要障礙[1]。重金屬通常具有持久的生物有效性并且在土壤中有較長(zhǎng)的滯留時(shí)間，它不僅會(huì)積累在農(nóng)業(yè)土壤中，導(dǎo)致土壤肥力降低、土壤生物學(xué)退化和作物生產(chǎn)力降低[2]，而且也容易通過(guò)呼吸、皮膚接觸或攝入等暴露途徑進(jìn)入動(dòng)植物和人體內(nèi)，并通過(guò)食物鏈進(jìn)行生物放大，給食品安全和人類健康帶來(lái)巨大風(fēng)險(xiǎn)[3]。

目前，對(duì)土壤重金屬的研究主要關(guān)注其污染特征、空間分布和來(lái)源解析。其中重金屬空間分布研究通過(guò)地統(tǒng)計(jì)方法來(lái)預(yù)測(cè)未采樣位置的值。地統(tǒng)計(jì)學(xué)目前被廣泛應(yīng)用于土壤、環(huán)境和生態(tài)學(xué)的研究當(dāng)中，它是識(shí)別重金屬污染區(qū)域的有效、直觀方法[4-5]?？死锝穑↘riging）法是目前地統(tǒng)計(jì)學(xué)中應(yīng)用最廣的方法之一，其原理是通過(guò)計(jì)算預(yù)測(cè)點(diǎn)與采樣點(diǎn)的空間關(guān)聯(lián)性，并使估計(jì)方差最小化來(lái)對(duì)預(yù)測(cè)位置進(jìn)行無(wú)偏估計(jì)[1]。然而，大多數(shù)地統(tǒng)計(jì)方法（包括克里金法）并未完全考慮對(duì)空間預(yù)測(cè)很重要的先驗(yàn)信息，并且難以輕易處理能夠提供輔助信息的大量相關(guān)知識(shí)（如同時(shí)處理數(shù)值變量和分類變量），因而限制了它們組合各種數(shù)據(jù)進(jìn)行更準(zhǔn)確的估計(jì)的能力[6]。

本文采用貝葉斯最大熵（Bayesian maximum entropy,BME）[7]將土壤母質(zhì)類型作為輔助信息，來(lái)提高土壤重金屬的空間預(yù)測(cè)精度。BME作為一種相對(duì)新穎的地統(tǒng)計(jì)學(xué)方法，能夠有效地綜合利用各種類型數(shù)據(jù)，從而提高預(yù)測(cè)的準(zhǔn)確性。這些數(shù)據(jù)能以不同的形式表達(dá)，如區(qū)間、概率密度函數(shù)（probability density function,PDF）等[8]。不同于其他傳統(tǒng)的統(tǒng)計(jì)方法，BME還具有以下優(yōu)點(diǎn)：1）不需要數(shù)據(jù)滿足基本的概率密度假設(shè)，即非高斯定律也可以應(yīng)用于統(tǒng)計(jì)分析；2）可以處理空間和時(shí)間非平穩(wěn)數(shù)據(jù)的預(yù)測(cè)，并可將物理規(guī)律和各種輔助信息（BME被稱之為軟數(shù)據(jù)）引入預(yù)測(cè)過(guò)程；3）能得到不同的分布圖（如期望值、95%預(yù)測(cè)區(qū)間、殘差等）[9-10]。本文在描述研究區(qū)土壤重金屬基本污染特征的基礎(chǔ)上，重點(diǎn)對(duì)比了BME與普通克里金（ordinary Kriging,OK）法對(duì)土壤重金屬的預(yù)測(cè)精度，為綜合利用輔助信息來(lái)提高土壤重金屬預(yù)測(cè)精度提供新的研究思路。

1 材料與方法

1.1 研究區(qū)概況

本文選取浙江省杭州市為研究區(qū)。杭州位于中國(guó)東南沿海（29°11′～30°33′N，118°21′～120°30′E），是浙江省的省會(huì)城市，也是長(zhǎng)三角城市群中的重要城市之一，面積約16 596 km2，總?cè)丝诩s919萬(wàn)。杭州屬于亞熱帶季風(fēng)氣候區(qū)，夏季炎熱潮濕，冬季寒冷干燥；年平均氣溫和相對(duì)濕度分別為17.8℃和70.3%，年均降水量和日照時(shí)數(shù)分別為1 454 mm和1 765 h。其東北部屬浙江北部平原，海拔低平，土壤肥沃，河網(wǎng)密布，是杭嘉湖平原和蕭紹平原的重要組成部分；西南部屬浙西中低山丘陵，地形復(fù)雜。

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

根據(jù)當(dāng)?shù)夭煌寥滥纲|(zhì)的分布和面積情況，按面積等因素確定不同母質(zhì)類型所需要的采樣點(diǎn)數(shù)，同時(shí)結(jié)合第二次土壤普查的相關(guān)資料及實(shí)地調(diào)查，綜合考慮海拔、地形地貌、土地利用類型及采樣點(diǎn)的代表性和可操作性等因素，預(yù)布設(shè)500個(gè)采樣點(diǎn)。在實(shí)地采樣過(guò)程中，根據(jù)預(yù)設(shè)點(diǎn)周邊環(huán)境進(jìn)行適當(dāng)調(diào)整，采用全球定位系統(tǒng)（global position system,GPS）記錄采樣點(diǎn)實(shí)際位置，形成采樣點(diǎn)分布圖（圖1）。每個(gè)采樣點(diǎn)的采樣深度為0～20 cm，在附近10 m范圍內(nèi)采集5個(gè)子樣本，充分混合，從而獲得具有代表性的土壤樣本。將土壤樣品在室溫下風(fēng)干，研磨至過(guò)100目網(wǎng)篩，進(jìn)行重金屬元素全量化學(xué)分析。土壤中重金屬含量采用標(biāo)準(zhǔn)方法測(cè)定，其中：鉻（Cr）、鉛（Pb）、鎘（Cd）和砷（As）通過(guò)硝酸-高氯酸混合消解，采用電感耦合等離子體質(zhì)譜（inductively coupled plasma mass spectrometry,ICP-MS）法測(cè)定[11]；汞（Hg）通過(guò)硝酸-過(guò)氧化氫微波混合消解，采用原子熒光光譜法測(cè)定[12]。通過(guò)加入隨機(jī)重復(fù)樣品和標(biāo)準(zhǔn)樣品對(duì)重金屬含量測(cè)定過(guò)程進(jìn)行質(zhì)量控制[13]，標(biāo)準(zhǔn)樣品的回收率在90%～110%之間，重復(fù)樣品的相對(duì)標(biāo)準(zhǔn)偏差在3%～8%之間。

杭州市土壤母質(zhì)二級(jí)分類數(shù)據(jù)提取自浙江省土壤數(shù)據(jù)庫(kù)[14]，該數(shù)據(jù)庫(kù)是以浙江省第二次土壤普查中編制的1∶5萬(wàn)比例尺土壤圖件為基礎(chǔ)，在搜集、整理數(shù)據(jù)資料的基礎(chǔ)上，進(jìn)行圖形掃描、矢量化處理，按照統(tǒng)一的數(shù)據(jù)庫(kù)建設(shè)標(biāo)準(zhǔn)而建成的具有正確拓?fù)潢P(guān)系的土壤數(shù)據(jù)庫(kù)，其分布如圖1所示。

圖1 杭州市土壤母質(zhì)及重金屬采樣點(diǎn)分布圖Fig.1 Distribution of soil parent materials and sampling points of heavy metals in Hangzhou City

1.3 空間預(yù)測(cè)方法

1.3.1 普通克里金插值法

作為一種傳統(tǒng)的地統(tǒng)計(jì)預(yù)測(cè)技術(shù)，OK已成功應(yīng)用于土壤科學(xué)，特別是對(duì)土壤重金屬的空間預(yù)測(cè)中[1,11]。它是在有限的離散數(shù)據(jù)基礎(chǔ)上實(shí)現(xiàn)對(duì)空間數(shù)據(jù)最優(yōu)線性無(wú)偏估計(jì)的工具，是探索變量在空間上的分布結(jié)構(gòu)和變異特性的有效研究方法[15]。它通過(guò)半方差函數(shù)來(lái)計(jì)算并模擬采樣點(diǎn)的可變性，從而預(yù)測(cè)未知位置的值，并可通過(guò)多種模型（高斯、指數(shù)、球形等）進(jìn)行擬合[16]。經(jīng)驗(yàn)半方差函數(shù)計(jì)算公式如下：

1.3.2 貝葉斯最大熵法

BME在時(shí)空隨機(jī)場(chǎng)的理論框架[17]下同時(shí)考慮多種數(shù)據(jù)的自相關(guān)和相互關(guān)系，支持多種數(shù)據(jù)表示方式（如區(qū)間、PDF等）來(lái)引入軟數(shù)據(jù)（輔助信息），從而提高統(tǒng)計(jì)分析和預(yù)測(cè)制圖的精度[7-9]。BME主要分為先驗(yàn)、后驗(yàn)和成圖3個(gè)階段。

1）先驗(yàn)階段：基于現(xiàn)有的一般知識(shí)（如自然定律、數(shù)學(xué)期望、半方差等），通過(guò)最大熵原理最大化期望的信息數(shù)據(jù)。熵通常用于描述系統(tǒng)的不確定性，公式為：

式中：S(X)和S(p1,p2,…,pn)是熵的函數(shù)；pi是事件i發(fā)生的概率。當(dāng)引入連續(xù)變量m(x)時(shí)，熵的函數(shù)可以變形為：

2）后驗(yàn)階段：綜合研究區(qū)的具體數(shù)據(jù)并基于貝葉斯條件概率原理，對(duì)第一步中的先驗(yàn)概率密度函數(shù)進(jìn)行更新，以獲得預(yù)測(cè)位置的后驗(yàn)概率密度函數(shù)。這一階段應(yīng)用的數(shù)據(jù)可以分為兩大類：硬數(shù)據(jù)和軟數(shù)據(jù)。本研究中的硬數(shù)據(jù)是指通過(guò)化學(xué)分析得到的土壤重金屬含量，具有高置信度，軟數(shù)據(jù)是指與土壤重金屬含量有關(guān)的土壤母質(zhì)類型信息[18-19]。后驗(yàn)概率密度函數(shù)計(jì)算公式如下：

式中：f(Ix0)是預(yù)測(cè)位置x0處的概率密度函數(shù)；Ihard和Rsoft分別是預(yù)測(cè)點(diǎn)一定范圍內(nèi)的硬數(shù)據(jù)和軟數(shù)據(jù)。在本研究中土壤母質(zhì)類型以區(qū)間方式的軟數(shù)據(jù)引入，所以后驗(yàn)概率密度函數(shù)可以轉(zhuǎn)換為：

土壤母質(zhì)類型具體的引入方式為：將采樣點(diǎn)根據(jù)土壤重金屬含量M等分為n組，將土壤母質(zhì)根據(jù)不同的類型劃分為不同的種類Z，記錄預(yù)測(cè)點(diǎn)x0處對(duì)應(yīng)的土壤母質(zhì)所屬的類型Zi，那么x0處根據(jù)土壤母質(zhì)類型所計(jì)算的軟數(shù)據(jù)概率為

式中：Rsoft是軟數(shù)據(jù)的概率函數(shù)；Counti是屬于土壤母質(zhì)類型Zi的樣點(diǎn)個(gè)數(shù)；Count(n)i是同時(shí)屬于土壤母質(zhì)類型Zi和重金屬含量分組Mn的樣點(diǎn)個(gè)數(shù)。

3）成圖階段：根據(jù)以上2步計(jì)算，由于得到的預(yù)測(cè)結(jié)果是概率密度函數(shù)，所以根據(jù)研究目的可以生成不同的空間分布圖（如期望、眾數(shù)、中位數(shù)、95%置信區(qū)間等）。

1.4 預(yù)測(cè)結(jié)果比較指標(biāo)

通過(guò)十折交叉驗(yàn)證（將數(shù)據(jù)集等分成10份，輪流將其中9份作為訓(xùn)練集，1份作為驗(yàn)證集，10次試驗(yàn)結(jié)果的正確率的平均值作為對(duì)算法精度的估計(jì)）檢驗(yàn)OK和BME插值方法對(duì)土壤重金屬空間預(yù)測(cè)的精度。以平均絕對(duì)誤差（mean absolute error,MAE）和均方根誤差（root mean square error,RMSE）2個(gè)指標(biāo)來(lái)進(jìn)行比較，計(jì)算公式如下：

式中：n為樣點(diǎn)的個(gè)數(shù)；yi和xi分別表示預(yù)測(cè)值和實(shí)際值。MAE和RMSE越小，表明該方法的預(yù)測(cè)精度越高。

2 結(jié)果與討論

2.1 土壤重金屬含量描述性特征統(tǒng)計(jì)

表1展示了土壤中重金屬含量的基本統(tǒng)計(jì)特征及浙江省土壤環(huán)境背景值[20]和國(guó)家《土壤環(huán)境質(zhì)量標(biāo)準(zhǔn)》（GB 15618—1995）的二級(jí)標(biāo)準(zhǔn)值。As、Cd、Cr、Hg和Pb的質(zhì)量分?jǐn)?shù)范圍為2.63～43.20、0.05～1.41、10.32～104.00、0.03～0.50和16.90～62.60 mg/kg，其平均值分別為9.00、0.27、52.88、0.13和31.67 mg/kg，均低于二級(jí)土壤環(huán)境質(zhì)量標(biāo)準(zhǔn)值，表明重金屬不存在超標(biāo)現(xiàn)象。但與浙江當(dāng)?shù)丨h(huán)境背景值相比，盡管Cr、Hg和Pb的平均含量低于相應(yīng)的背景值，但是As和Cd的平均含量都高于背景值，分別是背景值的1.31倍和1.59倍，這可能是由于過(guò)去幾十年的工農(nóng)業(yè)生產(chǎn)活動(dòng)和城市化進(jìn)程所造成的污染積累[18,20]。各重金屬的變異系數(shù)排序?yàn)镃d（80%）＞As（69%）＞Hg（61%）＞Cr（38%）＞Pb（28%），Cd、As和Hg較高的變異系數(shù)表明其在空間分布上存在很大的差異性，說(shuō)明它們受到了人類活動(dòng)的強(qiáng)烈干擾。此外，根據(jù)偏度和峰度值及Kolmogorov-Smirnov正態(tài)性檢驗(yàn)（P＜0.01），各重金屬含量均呈非正態(tài)分布，所以在空間插值前對(duì)重金屬含量進(jìn)行了標(biāo)準(zhǔn)正態(tài)轉(zhuǎn)換。

表1 杭州市土壤重金屬含量描述性統(tǒng)計(jì)Table 1 Descriptive statistics of soil heavy metals in Hangzhou City

2.2 半方差模型比較

圖2展示了用OK和BME法對(duì)不同重金屬計(jì)算得到的經(jīng)驗(yàn)半方差（圓圈）和擬合得到的理論半方差模型（線條）。對(duì)比OK和BME法的半方差模型發(fā)現(xiàn)：用2種計(jì)算方法得到的Cr的塊金效應(yīng)結(jié)果相同，均為0，表明Cr存在很強(qiáng)的空間關(guān)聯(lián)性，隨機(jī)誤差較小，受人為因素干擾較小。其余4種重金屬在加入土壤母質(zhì)信息的BME模型中其塊金效應(yīng)均小于OK模型，表明BME模型對(duì)重金屬的空間結(jié)構(gòu)具有更強(qiáng)的解釋能力。此外，從OK模型中可以看出，Cd、As、Hg和Pb的塊金效應(yīng)較大，表明這類重金屬存在較大的隨機(jī)誤差，其來(lái)源受人為因素干擾較大。

2.3 預(yù)測(cè)結(jié)果對(duì)比

由OK和BME法產(chǎn)生的重金屬預(yù)測(cè)殘差分布如圖3所示?？傮w上，BME對(duì)于重金屬空間預(yù)測(cè)的精度高于OK，表現(xiàn)為在0值處有更高的頻率分布和更窄的殘差區(qū)間。這種結(jié)果也得到表2中交叉驗(yàn)證下的MAE和RMSE證實(shí)。相較于OK，BME對(duì)Cr、Pb、Cd、Hg和As預(yù)測(cè)的MAE分別降低了13.82%、15.55%、1.87%、8.25%和19.45%，同時(shí)，RMSE分別降低了14.34%、19.72%、6.22%、9.56%和39.42%。由以上對(duì)比可知，與OK方法比較，利用土壤母質(zhì)類型作為輔助變量的BME方法能有效提高土壤重金屬預(yù)測(cè)精度。

圖2 克里金和貝葉斯最大熵半方差模型對(duì)比Fig.2 Semi-variance model comparison of Kriging and Bayesian maximum entropy methods

圖3 貝葉斯最大熵（BME）和克里金（OK）法預(yù)測(cè)的殘差分布對(duì)比圖Fig.3 Predicted residual distribution comparison of Bayesian maximum entropy(BME)and Kriging(OK)methods

表2 土壤重金屬BME和OK預(yù)測(cè)方法的計(jì)算精度對(duì)比Table 2 Comparison of estimation accuracy of Bayesian maximum entropy(BME)and Kriging(OK)methods

2.4 土壤重金屬空間分布

用BME法預(yù)測(cè)產(chǎn)生的土壤重金屬空間分布見(jiàn)圖4。從中可知，Pb和Hg具有類似的空間分布模式：均在東北部的城市地區(qū)重金屬含量較高?？焖俪鞘谢凸I(yè)化過(guò)程中造成的污染排放及交通運(yùn)輸帶來(lái)的污染物沉降可能是當(dāng)?shù)豍b和Hg含量較高的原因[21-23]。Cd和As也表現(xiàn)出相類似的空間分布規(guī)律：在西南部和中西部農(nóng)村地區(qū)重金屬含量較高。與化肥和農(nóng)藥密集使用密切相關(guān)的Cd，通常被視為農(nóng)業(yè)活動(dòng)的標(biāo)志性元素[24]，所以農(nóng)村地區(qū)的農(nóng)業(yè)活動(dòng)可能是造成當(dāng)?shù)刂亟饘俸可叩脑騕12,25]。此外，不同于As，Cd在中東部還存在一塊明顯的高含量區(qū)域，這與當(dāng)?shù)氐V業(yè)活動(dòng)密切相關(guān)[26]。由前文的半方差模型可知，Cr具有極好的空間結(jié)構(gòu)性，同時(shí)其變異系數(shù)也相對(duì)較低，說(shuō)明Cr的來(lái)源更主要是由于自然作用，受人為因素干擾較少。

3 結(jié)論與展望

本研究以杭州市土壤重金屬含量為研究對(duì)象，將OK和結(jié)合土壤母質(zhì)信息的BME方法進(jìn)行對(duì)比分析，比較其對(duì)重金屬含量預(yù)測(cè)的準(zhǔn)確性。由于BME法結(jié)合了土壤母質(zhì)輔助信息，它能夠更好地?cái)M合重金屬的空間變異規(guī)律，對(duì)重金屬預(yù)測(cè)的精度高于OK，具有更優(yōu)的MAE、RMSE值和殘差頻率分布。Pb和Hg在杭州東北部的城市地區(qū)具有較高的分布?？焖俪鞘谢凸I(yè)化過(guò)程中造成的污染排放及交通運(yùn)輸帶來(lái)的污染物沉降可能是其污染來(lái)源。Cd和As在西南部和中西部農(nóng)村地區(qū)存在明顯較高的分布，農(nóng)業(yè)活動(dòng)可能是其污染源。Cd在中東部還存在含量較高的區(qū)域，這與當(dāng)?shù)氐V業(yè)活動(dòng)密切相關(guān)。

圖4 基于貝葉斯最大熵（BME）預(yù)測(cè)的杭州市土壤重金屬空間分布圖Fig.4 Spatial distribution of soil heavy metals based on Bayesian maximum entropy(BME)prediction in Hangzhou City

本研究?jī)H采用了土壤母質(zhì)作為輔助信息來(lái)提高預(yù)測(cè)精度。考慮到BME法的多源數(shù)據(jù)融合能力，在今后的研究中還可以將與土壤重金屬密切相關(guān)的其他土壤屬性信息（如有機(jī)質(zhì)、pH等）、人為因素信息（工業(yè)、農(nóng)業(yè)、交通等數(shù)據(jù)）作為輔助數(shù)據(jù)引入，進(jìn)一步提高預(yù)測(cè)能力。同時(shí)，考慮到能夠使用輔助數(shù)據(jù)的克里金方法（如回歸克里金）對(duì)于土壤母質(zhì)這類面狀分類變量的引入存在一定的困難，本文只比較了BME與OK的差異，在后續(xù)引入數(shù)值變量的研究中可以進(jìn)一步比較其他克里金方法和BME法對(duì)于輔助數(shù)據(jù)的利用效率。