基于數(shù)學(xué)建模的城市表層土壤重金屬污染淺析

何雨森，李驊錦，蘆 慧，王玉蘭，魏友華

（成都理工大學(xué)，成都 610059）

0 前言

近年來(lái)，城市人口不斷加劇，城市經(jīng)濟(jì)快速發(fā)展，隨之產(chǎn)生的便是日益加劇的土壤地質(zhì)環(huán)境污染問(wèn)題，而導(dǎo)致土壤污染的原因主要?dú)w結(jié)于地質(zhì)活動(dòng)和人類(lèi)活動(dòng)。在現(xiàn)代相當(dāng)多的關(guān)于重金屬表層土壤污染物的前沿研究當(dāng)中，運(yùn)用數(shù)學(xué)模型（如對(duì)流模型，熱傳導(dǎo)方程模型等）針對(duì)動(dòng)態(tài)擴(kuò)散污染物，并基于GIS等計(jì)算機(jī)軟件模擬傳播過(guò)程并使之可視化的技術(shù)日趨成熟。作者在本文對(duì)于城市表層土壤重金屬污染問(wèn)題的探究，主要根據(jù)表層土壤（0cm～10cm）的采樣數(shù)據(jù)，重點(diǎn)分析解決人類(lèi)活動(dòng)對(duì)城市環(huán)境質(zhì)量的突出影響。

1 基于克里格插值的重金屬污染空間分析

克里格方法［1］（Kriging）公式可表示為式（1）。

式中Z（x0）為未知樣點(diǎn)的值；Z（xj）為未知樣點(diǎn)周?chē)阎獦颖军c(diǎn)的值；ωj為第j個(gè)已知樣本點(diǎn)對(duì)未知樣本點(diǎn)的權(quán)重；n為已知樣本點(diǎn)的個(gè)數(shù)。同時(shí)，作者引入了變差函數(shù)，來(lái)對(duì)目標(biāo)函數(shù)賦予權(quán)重：

由此，根據(jù)該模型進(jìn)行插值，可擬合出重金屬元素的空間分布圖。其流程圖見(jiàn)下頁(yè)圖1。

首先，將數(shù)據(jù)單位進(jìn)行統(tǒng)一，均以微克（μg）作為測(cè)量單位。隨后，利用克里格插值法加密數(shù)據(jù)點(diǎn)，以得到八種重金屬的空間分布圖，見(jiàn)下頁(yè)圖2。在圖2中，橫坐標(biāo)、縱坐標(biāo)為城市中各個(gè)功能區(qū)的位置，豎坐標(biāo)為各個(gè)元素在該位置的含量。

2 評(píng)價(jià)模型建立

2.1 單因子污染指數(shù)法建立評(píng)價(jià)指數(shù)模型

參照重金屬元素的背景值以及《中國(guó)土壤環(huán)境質(zhì)量標(biāo)準(zhǔn)》（GB15618－1995）［2］，得到了各種金屬的分級(jí)臨界值：

Xa（一級(jí)標(biāo)準(zhǔn)的最高值）。

Xc（二級(jí)標(biāo)準(zhǔn)的最高值）。

Xp（三級(jí)標(biāo)準(zhǔn)的最高值）。

作者采用單因子評(píng)價(jià)法［3］，建立重金屬污染評(píng)價(jià)指數(shù)模型見(jiàn)公式（2）。

對(duì)比標(biāo)準(zhǔn)得到了各個(gè)區(qū)域的各種重金屬的污染指數(shù)評(píng)級(jí)。結(jié)果如表1所示（見(jiàn)下頁(yè)）。

由于部份重金屬的含量存在極端情況，為了更嚴(yán)格地確定重金屬污染程度，作者在此采用單指標(biāo)最差原則［4］，工業(yè)區(qū)的Hg、Cu指標(biāo)已接近0.5，交通區(qū)的Cd、Cu、Zn 指標(biāo)也已接近0.5，占絕大比重。因此，針對(duì)工業(yè)區(qū)和交通區(qū)，我們修正的結(jié)果如表2所示（見(jiàn)下頁(yè)）。

表1 各區(qū)域的各種重金屬的污染指數(shù)評(píng)級(jí)Tab.1 The regional various heavy metal pollution index rating

表2 修訂后該城區(qū)不同區(qū)域重金屬污染級(jí)別Tab.2 Revised heavy metal pollution level in different region

2.2 各重金屬元素含量對(duì)區(qū)域污染程度的隸屬度評(píng)價(jià)

為了得到該城區(qū)中各個(gè)重金屬污染元素對(duì)土壤的綜合影響，對(duì)于每一個(gè)重金屬元素，需分別計(jì)算出每個(gè)功能區(qū)的相應(yīng)因素，對(duì)于各個(gè)等級(jí)的隸屬度，其步驟大致如下：

（1）建立因子集U＝｛As、Cd、Cr、Cu、Hg、Ni、Pb、Zn｝。

（2）計(jì)算每種金屬的隸屬度函數(shù)為式（3）。

通過(guò)分析表3中的數(shù)據(jù)可以得出：

（1）生活區(qū)主要的污染因素為Pb、Cu。

（2）工業(yè)區(qū)主要污染因素為Cu、Hg、Pb、Zn。

（3）山區(qū)均符合標(biāo)準(zhǔn)。

（4）交通區(qū)主要污染因素為Cr、Pb。

（5）公園綠地區(qū)主要污染因素為Cd、Pb、Cu、Hg。

2.3 因子的灰色關(guān)聯(lián)分析

如果重金屬之間存在顯著關(guān)聯(lián)，則其來(lái)源可能相同，否則來(lái)源可能不止一個(gè)。為此，作者建立了關(guān)聯(lián)度模型［5］，研究各個(gè)污染物因子之間的關(guān)聯(lián)性。

令ZBi為第i個(gè)重金屬元素，則有：ZBi、ZB0在第k類(lèi)地區(qū)的關(guān)聯(lián)系數(shù)，如公式（4）。

其中 △i（k）＝｜xki－xoi｜，（i＝1，2，…，m；k＝1，2，…，n）；a＝min△i（k）；b＝max。

ZBi與ZB0的總體關(guān)聯(lián)度為：

從而得到不同功能區(qū)各個(gè)重金屬元素的關(guān)聯(lián)度系數(shù)。作者以Zn在各個(gè)區(qū)域和其它金屬元素的關(guān)聯(lián)度為例，結(jié)果見(jiàn)下頁(yè)表4。

經(jīng)過(guò)關(guān)聯(lián)度分析，最終得出該城區(qū)各功能區(qū)重金屬污染均與工業(yè)污染有關(guān)，且各個(gè)重金屬元素之間有很強(qiáng)的關(guān)聯(lián)性。因此，該城區(qū)工業(yè)污染相對(duì)嚴(yán)重。

3 模擬污染物的傳播及確定污染物的位置

3.1 基于熱傳導(dǎo)方程的傳播過(guò)程模擬

經(jīng)查詢(xún)文獻(xiàn)發(fā)現(xiàn)，在平面上重金屬污染物的傳播和平面上熱傳導(dǎo)過(guò)程有相似性［6］，故利用熱傳導(dǎo)過(guò)程的方法與思路模擬其傳播過(guò)程。

2.2 兩組患者治療前后胎兒生長(zhǎng)指標(biāo)比較 治療前，兩組患者孕20周時(shí)胎兒BPD及FL比較，差異無(wú)統(tǒng)計(jì)學(xué)意義(P>0.05)；治療后，用藥8周后胎兒BPD及FL顯著高于治療前(P<0.05)，且治療組明顯高于對(duì)照組(P<0.05)。見(jiàn)表2。

首先假設(shè)在某一個(gè)小區(qū)域內(nèi)僅有一個(gè)污染點(diǎn)“0”，在污染源處的濃度會(huì)很高，即為濃度的唯一極值點(diǎn)。此點(diǎn)類(lèi)似于熱傳導(dǎo)模型中的點(diǎn)熱源，周?chē)鷧^(qū)域的含量大小也類(lèi)似于點(diǎn)熱源周?chē)氖軣岢潭取?/p>

傳統(tǒng)的熱傳導(dǎo)方程公式［7］如公式（5）。

3.2 差分法網(wǎng)格化

熱傳導(dǎo)方程是典型的偏微分方程，通過(guò)偏微分方程數(shù)值解中的差分法，可以設(shè)計(jì)疊代算法程序計(jì)算每種重金屬在城區(qū)的分布。首先，根據(jù)實(shí)際情況設(shè)定邊界、初始條件；然后，離散化方程；最后，在程序上計(jì)算得到一系列在設(shè)定的坐標(biāo)系下網(wǎng)格結(jié)點(diǎn)處的值，即為重金屬元素的濃度值。

表3 八種重金屬元素對(duì)于五個(gè)功能區(qū)污染狀況隸屬度Tab.3 Degree of membership of 8kinds of heavy metals on 5functional areas pollution

表4 Zn在各個(gè)功能區(qū)與其它金屬的關(guān)聯(lián)度Tab.4 Association between Zn and other metal element in various functional areas

3.3 模擬退火算法找到污染源

通過(guò)差分法，已經(jīng)基本得到了各種元素在城區(qū)空間的靜態(tài)分布狀況，接下來(lái)就是尋找污染源。作者運(yùn)用模擬退火算法，通過(guò)構(gòu)建模型尋找各種重金屬元素濃度的最高值，來(lái)確定污染源的位置。

在使用退火算法的過(guò)程當(dāng)中，默認(rèn)每次迭代得到的極大值即為最優(yōu)解。

退火系統(tǒng)在早期允許差解存在，存在概率為P，這個(gè)概率會(huì)越來(lái)越小。由于目標(biāo)是要找到最優(yōu)解，所以系統(tǒng)后期基本不接受差解，可以用公式（6）來(lái)表示。

為了尋找最優(yōu)解，對(duì)每種重金屬元素的數(shù)據(jù)進(jìn)行十次測(cè)試，取十次運(yùn)算的最佳值作為這組數(shù)據(jù)的實(shí)驗(yàn)結(jié)果。利用Matlab軟件，按照上述算法編程實(shí)現(xiàn)，最終得到每種元素在該城區(qū)內(nèi)極大值點(diǎn)的坐標(biāo)，見(jiàn)后面表5。

3.4 結(jié)果驗(yàn)證分析

為了提高模型的科學(xué)性、準(zhǔn)確性，需要對(duì)模擬結(jié)果進(jìn)行驗(yàn)證。作者將各種金屬的模擬來(lái)源位置與城市區(qū)域圖進(jìn)行疊加，如圖4所示（見(jiàn)下頁(yè)）。

通過(guò)觀察表5可知，As與Cu、Cr與Ni、Zn來(lái)源于同一位置，且這八種元素的位置距離較近。結(jié)合數(shù)據(jù)，極大值坐標(biāo)位于（6，8）的有As、Cd、Cu，它們的污染源位置主要在工業(yè)區(qū)；而極大值坐標(biāo)位于（9，22）的有Cr、Hg、Ni、Zn，它們的污染源位置主要在交通區(qū)；極大值坐標(biāo)位于（16）的即Pb，它的污染源位置在生活區(qū)。

結(jié)合關(guān)聯(lián)度模型，可知各種重金屬元素具有很強(qiáng)的關(guān)聯(lián)性，可能是由于同一污染源所致。并且，通過(guò)對(duì)關(guān)聯(lián)度模型中污染主要原因分析可知，該城區(qū)主要受工業(yè)污染的干擾導(dǎo)致土壤質(zhì)量不斷惡化。可以推斷，該城區(qū)污染源的位置在坐標(biāo)（1647＊895～6534＊6018）這一片區(qū)域內(nèi)，且主要由坐標(biāo)（1647，2728）和坐標(biāo)（2383，3692）這兩地的工業(yè)區(qū)所致。因此坐標(biāo)（1647，2728）和坐標(biāo)（2383，3692）這兩地為該城區(qū)污染源的位置。

圖3 模擬退火模型的基本算法Fig.3 Basic algorithms of simulated annealing module

圖4 各種金屬的模擬來(lái)源位置與城市區(qū)域疊加圖Fig.4 Layout chart combined with urban area and location of simulated source of various metal element

表5 每種重金屬元素在該城區(qū)內(nèi)極大值點(diǎn)坐標(biāo)Tab.5 Maximum point coordination of each kind of heavy metal element in the area

4 結(jié)論

綜上所述，作者所使用的方法以及所獲得的結(jié)論，可以為其它城市重金屬表層土壤污染相關(guān)的研究與分析提供參考。

［1］秦偉良，劉悅.空間插值法在降水分布中的應(yīng)用［J］.南京信息工程大學(xué)學(xué)報(bào)：自然科學(xué)版.2010（2）：162.

［2］百度百科.土壤環(huán)境質(zhì)量標(biāo)準(zhǔn).http：／／baike.baidu.com／view／1738501.htm.2011.

［3］王芬，彭國(guó)照，蔣錦剛.基于雙層組合神經(jīng)網(wǎng)絡(luò)的川芎道地產(chǎn)區(qū)土壤重金屬污染評(píng)價(jià)與GIS可視化［J］.成都信息工程學(xué)院資源環(huán)境學(xué)院學(xué)報(bào)，2010（3）：75.

［4］段雪梅，蔡煥興，巢文軍.南京市表層土壤重金屬污染特征及污染來(lái)源［J］.南京市環(huán)境科學(xué)與管理學(xué)報(bào)，2010（10）：31.

［5］劉光廷，邱德隆.三維熱傳導(dǎo)問(wèn)題的間接邊界元法［J］.清華大學(xué)學(xué)報(bào)：自然科學(xué)版，2010（36）：8.

［6］王學(xué)賀.離散正則方法在一維熱傳導(dǎo)方程尋源反問(wèn)題中的應(yīng)用［J］.燕山大學(xué)學(xué)報(bào)，2008（32）：88.

［7］謝紅燕.基于并行模擬退火算法的VRP 問(wèn)題研究［J］.物流技術(shù)月刊，2010（8）：67.

［8］湯世明.基于GIS的大氣污染物擴(kuò)散模擬分析系統(tǒng)設(shè)計(jì)與實(shí)現(xiàn)［J］.水電能源科學(xué)，2011（9）：29.

［9］陳寶.上海第4層淤泥質(zhì)粉質(zhì)粘土的熱傳導(dǎo)特性［J］.地下空間與工程學(xué)報(bào)，2011（8）：25.

［10］劉勇，岳玲玲，李晉昌.太陽(yáng)市土壤重金屬污染及其潛在生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)［J］.環(huán)境科學(xué)學(xué)報(bào)，2011，222（6）：67.

［11］唐發(fā)靜，祖艷群.土壤重金屬空間變異的研究方法［J］.云南農(nóng)業(yè)大學(xué)學(xué)報(bào)，2008（4）：76.

［12］LI X D，LEE S L，WONG S C，et al.The study of metal contamination in urban soils of Hong Kong using a GIS－based approcah［M］.Environmental Pollution.2004.

［13］HUANG KUO－MING，LIN SAULWOOD.Consequences and implication of heavy metal spatial variations in sediments of the Keelung River drainage basin，Taiwan［M］.Chemosphere.2003.

［14］王軍，陳振樓，王初，等.上海崇明島蔬菜地土壤重金屬含量與生態(tài)風(fēng)險(xiǎn)預(yù)警評(píng)估［J］.環(huán)境科學(xué)，2007（3）：133.