2011數(shù)學(xué)建模全國一等獎(jiǎng)A城市表層土壤重金屬污染分析

1、2011高教社杯全國大學(xué)生數(shù)學(xué)建模競賽承 諾 書我們仔細(xì)閱讀了中國大學(xué)生數(shù)學(xué)建模競賽的競賽規(guī)則.我們完全明白，在競賽開始后參賽隊(duì)員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊(duì)外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問題。我們知道，抄襲別人的成果是違反競賽規(guī)則的, 如果引用別人的成果或其他公開的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻(xiàn)的表述方式在正文引用處和參考文獻(xiàn)中明確列出。我們鄭重承諾，嚴(yán)格遵守競賽規(guī)則，以保證競賽的公正、公平性。如有違反競賽規(guī)則的行為，我們將受到嚴(yán)肅處理。賽區(qū)評閱編號（由賽區(qū)組委會評閱前進(jìn)行編號）：2011高教社杯全國大學(xué)生數(shù)學(xué)建模競賽編 號 專

2、用 頁賽區(qū)評閱編號（由賽區(qū)組委會評閱前進(jìn)行編號）：賽區(qū)評閱記錄（可供賽區(qū)評閱時(shí)使用）：評閱人評分備注全國統(tǒng)一編號（由賽區(qū)組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進(jìn)行編號）：城市表層土壤重金屬污染分析摘要本文主要研究重金屬對城市表層土壤污染的問題,我們根據(jù)題目所給定的一些數(shù)據(jù)和信息分析并建立了擴(kuò)散傳播模型、權(quán)重分配模型、對比模型和轉(zhuǎn)換模型解決問題。首先，我們利用Matlab軟件擬出該城區(qū)地勢圖（圖1），根據(jù)所給數(shù)據(jù)繪出該地區(qū)的三維地勢及采樣點(diǎn)在其上的綜合空間分布圖。之后將8種重金屬的濃度等高線投影到該地區(qū)三維地形圖曲面上，接著分別計(jì)算8種重金屬在五個(gè)區(qū)域的平均值，立體圖和平面圖（

3、圖1附件）相結(jié)合便可得出8種重金屬元素在該城區(qū)的空間分布。其次，在確定該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度時(shí)，我們運(yùn)用兩種方法進(jìn)行解答。先假設(shè)各重金屬毒性及其它性質(zhì)相同，運(yùn)用公式求出各區(qū)域各金屬相對于背景平均值的比值作為金屬污染程度，再運(yùn)用求出各區(qū)域重金屬污染程度，并將各區(qū)進(jìn)行比較。之后，我們加上各重金屬的毒性，對各重金屬求出權(quán)數(shù)，再結(jié)合國標(biāo)重金屬污染等級和已知的各組數(shù)據(jù)來確定金屬的污染程度。由上述兩種方法的對比，更準(zhǔn)確地得出重金屬對各區(qū)的影響程度。即： 工業(yè)區(qū)>交通區(qū)>生活區(qū)>公園綠地區(qū)>山區(qū)并根據(jù)第一個(gè)模型的數(shù)據(jù)來說明重金屬污染的主要原因。再次，對重金屬污染物的傳播特征

4、進(jìn)行了分析，判斷出重金屬污染物主要是通過大氣、土壤和水流進(jìn)行傳播。在分析之中，我們得出這三種狀態(tài)的傳播并不是孤立存在的，而是可以相互影響和疊加的，因此，我們分別建立三個(gè)傳播模型，再對這三個(gè)傳播模型進(jìn)行了時(shí)間和空間上的擬合，得出重金屬濃度最高的區(qū)域圖，并結(jié)合各重金屬的分布圖（圖6）來確定各污染源的位置。最后，本題中只給出了重金屬對土壤的污染，對于研究城市地質(zhì)環(huán)境的演變模式，還需要搜集一些信息（圖7）。根據(jù)每種因素對地質(zhì)環(huán)境的影響程度進(jìn)行由定性到定量的轉(zhuǎn)化。建立同一地質(zhì)時(shí)期地質(zhì)環(huán)境中各因素的正影響和負(fù)影響的權(quán)重分配模型，再對這些權(quán)重進(jìn)行驗(yàn)算和修正。從而，根據(jù)這些權(quán)重再建立預(yù)測模型便可反向推出各重金

5、屬對不同時(shí)期地質(zhì)環(huán)境的影響，得出隨時(shí)間變化的地質(zhì)環(huán)境的演變模式。結(jié)論：在本次模型建立中，我們得出以下結(jié)論：1.重金屬在各個(gè)區(qū)域中的污染嚴(yán)重程度為：工業(yè)區(qū)>交通區(qū)>生活區(qū)>公園綠地區(qū)>山區(qū)2.各重金屬的污染源主要分布在工業(yè)區(qū)和交通區(qū) 關(guān)鍵詞：重金屬污染 三維地形圖 時(shí)空結(jié)合 地質(zhì)環(huán)境演變 影響因子權(quán)重一問題重述1.問題背景目前，社會經(jīng)濟(jì)發(fā)展迅速，人口數(shù)量不斷增加，環(huán)境污染現(xiàn)象日顯突出，尤其重金屬對土壤的污染更受廣泛關(guān)注。土壤狀況直接影響著動植物的生長和安全，甚至通過食物鏈進(jìn)入人體，導(dǎo)致一些慢性疾病的發(fā)生。對于具有獨(dú)立的系統(tǒng)來說，人們的生活和生產(chǎn)將會給環(huán)境和土壤造成污染，而

6、且，每一個(gè)區(qū)域的功能不同，如山區(qū)、生活區(qū)、工業(yè)區(qū)、主干道路區(qū)和綠地區(qū)等，對環(huán)境和土壤的污染程度也不同。所以，做好調(diào)查分析，控制污染源是現(xiàn)今的關(guān)鍵。2.提出問題： （1）.根據(jù)題中所給各區(qū)域點(diǎn)的坐標(biāo)，繪制中該區(qū)的空間分布圖，計(jì)算不同重金屬對該區(qū)的影響。（2）.分析數(shù)據(jù)，考慮濃度分布和客觀因素，說明重金屬的來源。（3）.分析重金屬產(chǎn)生后的傳播特性，并建立模型確定污染源的位置。（4）.考慮重金屬濃度和客觀因素優(yōu)化模型并推廣，以便對城市地質(zhì)環(huán)境的演變做更準(zhǔn)確的分析和預(yù)測。二問題分析（1）.由于各種重金屬在各個(gè)區(qū)域的濃度都不同，所以對不同區(qū)域的影響程度也不同，根據(jù)題中所給采樣點(diǎn)的坐標(biāo)和各重金屬的濃度，用

7、MATLAB軟件繪制出重金屬在該區(qū)的空間分布圖及8種重金屬的濃度等高線在該地區(qū)三維地形圖曲面的投影圖。利用平均數(shù)法求得各種金屬在不同區(qū)域的影響。在考慮重金屬毒性的情況下，利用加權(quán)平均數(shù)求出各種金屬的污染程度。兩種方案作對比，找出最優(yōu)化方案。 （2）.在上題的結(jié)論基礎(chǔ)上，依據(jù)各種重金屬在不同區(qū)域的濃度大小，通過分析比較，并考慮各區(qū)域的特點(diǎn)，可以分析出重金屬的主要來源。 （3）.由于重金屬產(chǎn)生后的傳播途徑有很多，本題忽略其它因素，只考慮水體、大氣和土壤傳播，分別建立水體傳播模型，大氣沉降模型和固體傳播模型，并結(jié)合三者在傳播過程中在時(shí)間和空間上的相互影響和制約，對污染源的位置進(jìn)行準(zhǔn)確的定位。 （4）

8、.在研究城市地質(zhì)環(huán)境演變的過程中，除運(yùn)用到重金屬的濃度之外，還應(yīng)搜集重金屬的沉積和埋藏深度、生活污水，廢氣污染、地殼運(yùn)動植物和微生物的凈化。根據(jù)污染毒性的大小確定各項(xiàng)的權(quán)重，。由此建立模型推出不同時(shí)期地質(zhì)環(huán)境的變化，即演變模式。三問題假設(shè)1.假設(shè)重金屬元素在土壤和水中化學(xué)反應(yīng)均勻。2.假設(shè)各區(qū)域成土母質(zhì)中含重金屬的濃度是相同的。3.假設(shè)各地區(qū)重金屬分布穩(wěn)定，污染源排放量不變4.假設(shè)各數(shù)據(jù)擬合的函數(shù)是光滑曲線，且各外在因素都忽略。5.計(jì)算重金屬在空氣中傳播速度不考慮風(fēng)速過大的影響。四、 名詞解釋及符號說明1名詞解釋：1.地質(zhì)環(huán)境：包括巖石、水、氣和生物在內(nèi)的互相關(guān)聯(lián)的系統(tǒng)。 2.時(shí)空結(jié)合：將金屬

9、傳播的三種模型進(jìn)行擬合，找出污染源。2符號說明 25 i區(qū)域中j重金屬元素的濃度 i 區(qū)域中j重金屬元素的平均濃度 i 區(qū)域中調(diào)查重金屬的點(diǎn)數(shù) i 區(qū)域中j重金屬元素污染程度N 區(qū)域金屬對評價(jià)等級的隸屬度R 各金屬因子的權(quán)重構(gòu)成的向量D 各金屬因子對評價(jià)等級的隸屬度m 土壤環(huán)境質(zhì)量級別 金屬因子的隸屬度 某區(qū)域各金屬因子j在k級指標(biāo)P 各金屬因子的實(shí)測濃度 某地區(qū)第j個(gè)重金屬污染物的實(shí)測濃度; 第j個(gè)重金屬污染物的毒性級別指數(shù); 某地區(qū)第j個(gè)重金屬污染物的權(quán)重值 土顆粒表面所吸附的溶質(zhì)（j）的質(zhì)量分?jǐn)?shù) 孔隙中的溶質(zhì)的質(zhì)量濃度（）； 土壤孔隙比 土樣干密度（） 土壤孔隙率S 土壤容重 氣體密度

10、粒子密度r 粒子半徑 時(shí)間 阻滯因子） 擴(kuò)散系數(shù)，v為孔隙流體的流速 動力粘性系數(shù)W 粒子垂直方向速度 粒子受到的浮力 粒子受到的流體阻力 隨機(jī)力及大氣流動引起的隨機(jī)加速度五、 模型的建立及求解一、 重金屬元素的空間分布及污染程度模型1. 重金屬元素的空間分布：a、根據(jù)已知數(shù)據(jù)，我們運(yùn)用MATLAB 軟件將該地區(qū)的三維地勢及采樣點(diǎn)在其上的綜合空間分布圖繪制如下圖（1）所示在此圖中，我們能清晰的分出生活區(qū)、工業(yè)區(qū)、山區(qū)、交通區(qū)、公園綠地區(qū)。將生活區(qū)、工業(yè)區(qū)、山區(qū)、交通區(qū)、公園綠地區(qū)用i表示。i=(1、2、3、4、5) b、對八種重金屬元素的濃度進(jìn)行處理： 設(shè)八種重金屬元素的濃度用表示，As (g

11、/g)、Cd (ng/g)、Cr (g/g)、Cu (g/g)、Hg (ng/g)、 Ni (g/g)、 Pb (g/g)、 Zn (g/g)分別為j=(1、2、8) 對數(shù)據(jù)進(jìn)行處理分別求出各地區(qū)各金屬元素的平均濃度，根據(jù)公式 （1）計(jì)算得各數(shù)據(jù)如下表（1）所示： J I As (g/g)Cd (ng/g)Cr (g/g)Cu (g/g)Hg (ng/g)Ni (g/g)Pb (g/g)Zn (g/g)生活區(qū)6.27289.9669.0249.4093.0418.3469.11237.01工業(yè)區(qū)7.25393.1153.41127.54642.3619.8193.04277.93山區(qū)4.041

12、27.0027.5823.9930.0011.9357.4585.61交通區(qū)5.71360.0158.0562.21446.8217.6263.53242.85公園綠地區(qū)6.26280.5443.6430.19114.9915.2960.71154.24 c. 用Matlab軟件繪制8種重金屬的濃度等高線在該地區(qū)三維地形圖曲面的投影如下圖（圖1附件）： 結(jié)合圖（1）、圖（1）附件和表（1）清晰地看出8種主要重金屬元素在該城區(qū)各區(qū)域的空間分布。2. 重金屬元素污染程度模型：1)、在不考慮重金屬毒性的情況下：已知的背景區(qū)重金屬平均濃度如表（2）：8種主要重金屬元素的背景值元素平均值標(biāo)準(zhǔn)偏差范圍As

13、 (g/g)3.60.91.85.4Cd (ng/g)1303070190Cr (g/g)3191349Cu (g/g)13.23.66.020.4Hg (ng/g)3581951Ni (g/g)12.33.84.719.9Pb (g/g)3161943Zn (g/g)69144197設(shè)重金屬元素污染程度,則， (2)根據(jù)表（1）、表（2）的數(shù)據(jù)及公式(2)可得各區(qū)域的重金屬元素污染程度如表（3）jciAs (g/g)Cd (ng/g)Cr (g/g)Cu (g/g)Hg (ng/g)Ni (g/g)Pb (g/g)Zn (g/g)生活區(qū)1.742.362.233.742.661.492.23

14、3.43工業(yè)區(qū)2.013.021.729.6618.351.613.004.03山區(qū)1.120.980.891.820.860.971.851.24交通區(qū)1.592.771.874.7112.771.432.053.52公園綠地區(qū)1.742.161.412.893.291.241.962.24由于不考慮重金屬毒性，將各重金屬的污染性視為相同，則各地區(qū)金屬元素總體的污染程度即： （3）所以： 經(jīng)計(jì)算，可以得到以下關(guān)系式；>>>> ；即： 該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度工業(yè)區(qū)>交通區(qū)>生活區(qū)>公園綠地區(qū)>山區(qū)2)在考慮重金屬的毒性的情況下： 目前大

15、多數(shù)人在確定污染程度時(shí)僅僅考慮了重金屬污染物濃度超標(biāo)的情況,未考慮重金屬本身的毒性作用,這就有可能掩蓋有些濃度低但毒性大的有毒物的污染作用。 本模型法用隸屬度來描述模糊的污染分級界線,各評價(jià)等級的隸屬度再以各金屬因子的權(quán)重修正,則得到金屬因子對評價(jià)等級的隸屬度。則得到如下數(shù)學(xué)模型: （4） 式中：N 區(qū)域金屬對評價(jià)等級的隸屬度 R 各金屬因子的權(quán)重構(gòu)成的向量 D 各金屬因子對評價(jià)等級的隸屬度 a、 金屬因子隸屬度D的建立：為了運(yùn)算，我們需要建立隸屬度函數(shù)，并用隸屬度來描述土壤污染狀況的模糊界線。設(shè)土壤環(huán)境質(zhì)量分為m 個(gè)級別,則 V = (1,2, , m ) 在這里我們用降半梯形分布來刻畫隸屬

16、度: （5）式中： 金屬因子的隸屬度 某區(qū)域各金屬因子j在k級指標(biāo)(j = 1, 2, ,n; k = 1, 2, 3, m )中的標(biāo)準(zhǔn)值; P 各金屬因子的實(shí)測濃度 由此可得評價(jià)因子j對不同級別k的隸屬度矩陣D: (6)b、 金屬因子權(quán)重向量R 的確定: 將污染物濃度和毒性級別指數(shù)加權(quán)疊加,并作歸一化處理,得到某污染組分的權(quán)重公式: (7) (8)式中, 某地區(qū)第j個(gè)重金屬污染物的實(shí)測濃度; 第j個(gè)重金屬污染物的毒性級別指數(shù); 某地區(qū)第j個(gè)重金屬污染物的權(quán)重值,且 ，將各金屬因子的實(shí)測濃度值、毒性系數(shù)和選定的評價(jià)標(biāo)準(zhǔn)分別代入上式,可得到各金屬因子的權(quán)重值,由此組成某個(gè)區(qū)域各金屬因子的權(quán)重向量

17、:R = 、 , , , (9) 將權(quán)重向量R和隸屬度矩陣D帶入模型（4）可得區(qū)域金屬對評價(jià)等級的隸屬度NC、隸屬度函數(shù)的確定： 本文所選用的八種元素，選用國家土壤環(huán)境質(zhì)量標(biāo)準(zhǔn) 18 ( GB15618 -1995)進(jìn)行土壤環(huán)境評價(jià),土壤環(huán)境質(zhì)量標(biāo)準(zhǔn)見表（4）。根據(jù)表（1）和表（5）的數(shù)據(jù)，利用公式（5）計(jì)算各重金屬元素對應(yīng)于各土壤重金屬環(huán)境質(zhì)量等級的隸屬函數(shù),得到關(guān)系模糊矩陣。表（4）國家土壤環(huán)境質(zhì)量標(biāo)準(zhǔn)一級自然景值二級PH6.5 PH=6.57.5 PH 7.5三級PH6.5As (mg/kg)1530252030Cd (mg/kg)0.200.300.300.601.00Cr (mg/k

18、g)90150200250300Cu (mg/kg)3550100100400Hg (mg/kg)0.150.30.51.01.5Ni (mg/kg)40405060200Pb (mg/kg)35250300350500Zn (mg/kg)100200250300500 根據(jù)Hakanson制定的標(biāo)準(zhǔn)化重金屬毒性響應(yīng)系數(shù),分別對各重金屬對生物的毒性級別指數(shù)賦值：Zn=l， Cr=2，Cu=5，Pb=5，Ni=5， As=10，Cd=30，Hg=40，指數(shù)越大，代表毒性越大。按照（7）、（8）兩式，將表（1）的數(shù)據(jù)和毒性指數(shù)代入計(jì)算，得各區(qū)域各重金屬權(quán)重值，其中 R = 、 , , , = 0.

19、 0962 、0. 2838、 0. 0680、 0. 0585、 0. 3723、 0. 0211、0.0563、0.0538 （10）結(jié)合公式（4）、（6）、（10）及表（1）可得： =2.5158 =8.9537= 1.0546 =6.4033=2.5056所以：，即，工業(yè)區(qū)>交通區(qū)>生活區(qū)>公園綠地區(qū)>山區(qū)結(jié)合1）、2）兩個(gè)模型可以得到該城區(qū)各區(qū)域重金屬的污染程度為工業(yè)區(qū)>交通區(qū)>生活區(qū)>公園綠地區(qū)>山區(qū)二重金屬污染成因：根據(jù)上述兩個(gè)模型的結(jié)合與對比，我們可以得出下面的關(guān)系式：工業(yè)區(qū)>交通區(qū)>生活區(qū)>公園綠地區(qū)>山

20、區(qū)1.從整體的角度分析： 工業(yè)區(qū)污染最大，占整體區(qū)域的36%。分析其污染的主要來源有金屬礦山的開采、冶煉、重金屬尾礦、冶煉廢渣和礦渣堆放等。而且其擴(kuò)散和傳播污染的途徑也是非常廣泛的。 交通區(qū)重金屬的污染程度相對也較大，其污染主要體現(xiàn)為大氣污染，來自工業(yè)生產(chǎn)排放的廢氣、汽車尾氣排放生的大量含重金屬的有害氣體和粉塵。 生活區(qū)重金屬污染居中，重金屬的主要來源有生活垃圾、污水、醫(yī)療垃圾和細(xì)菌等。 公園綠地區(qū)重金屬污染程度較輕，其主要來源有農(nóng)業(yè)農(nóng)藥和化肥。但是，植物有吸收金屬礦物的作用，相對減輕了重金屬的危害。 山區(qū)重金屬污染最輕，只占到整體區(qū)域的8%，表明不僅重金屬的來源相對較少，而且，重金屬污染與重

21、力和海拔也有一定關(guān)系，山區(qū)中植物和一些微生物都有凈化的作用。2.從局部的角度分析： 從表（3）可以看出八種重金屬的主要分布區(qū)域，再根據(jù)這些區(qū)域的區(qū)別可以分析出各金屬污染物的成因。其主要來源如下表所示：污染物來源砷As礦山開采、冶煉、加工排放的三廢、飄塵、農(nóng)藥等鎘Cd電鍍、塑料電池電子、染料化工工業(yè)廢水鉻Cr電鍍、染料化工工業(yè)廢水、飄塵、礦山開采、冶煉、主要來源是汽車燃料和輪胎磨損。銅Cu電鍍工業(yè)廢水、農(nóng)藥汞Hg塑料電池電子工業(yè)廢水、礦山開采、冶煉、加工排放的三廢鎳Ni電鍍工業(yè)廢水、礦山開采、冶煉、加工排放的三廢鉛Pb電鍍工業(yè)廢水、礦山開采、冶煉、加工排放的三廢、汽車尾氣鋅Zn電鍍工業(yè)廢水三、重

22、金屬污染物傳播特性及污染源位置確定模型： 影響重金屬污染物傳播的因素有很多，如：地質(zhì)條件（如海拔，地殼運(yùn)動，土質(zhì)情況等）、氣候條件（如降雨量，氣溫，季風(fēng)等）、工業(yè)發(fā)展情況、人口密度、河流流向等。在這里，我們由于條件有限我們僅考慮地勢、河流等。重金屬污染物的傳播有許多種，但主要有三種方式,即：水體傳播、土壤傳播、大氣沉降。在此，我們將模型理想化，假設(shè)只有這三種方式，分別對其建立模型，之后再將三種模型進(jìn)行時(shí)間和空間上的幾何疊加。1.水體傳播模型的建立及求解：在此模型中，我們將問題理想化，認(rèn)為水體傳播僅僅是靠河流傳播的。假設(shè)河流沒有多大的改變。在河流中，重金屬污染物部分以可溶狀態(tài)存在，隨水流運(yùn)動，部

23、分以沙粒為載體，隨沙運(yùn)動，兩者不斷地發(fā)生化學(xué)和生物化學(xué)反應(yīng)，一般情況下，泥沙吸附解吸是重金屬污染物從不飽和液相轉(zhuǎn)移到固相的主要途徑。因此，在考慮重金屬在水環(huán)境中的物理遷移過程下，我們用水力學(xué)和泥沙運(yùn)動學(xué)的原理處理。假設(shè)，泥沙濃度和水相金屬離子濃度反應(yīng)分布均勻。在天然河流中, 水環(huán)境化學(xué)條件對某一河段、某一時(shí)段基本不變,變化較大并難處理的是水流泥沙因素。泥沙顆粒在水溶液中對重金屬污染物的吸附遵循蘭格茂爾吸附運(yùn)動學(xué)方程，將其對時(shí)間的偏微分 ，即對均勻沙有： （11）水流泥沙運(yùn)動數(shù)學(xué)模擬研究已有較長的歷史，在水流泥沙運(yùn)動確定后, 根據(jù)建立的重金屬遷移轉(zhuǎn)化方程, 吸附動力學(xué)方程及邊界、初始條件就可確定

24、河流重金屬污染物遷移轉(zhuǎn)化過程, 構(gòu)成河流一維重金屬污染物遷移轉(zhuǎn)化數(shù)學(xué)模型的基本方程如下：水流連續(xù)方程 （12） 水流運(yùn)動方程 （13）泥沙連續(xù)方程 （14）河床變形方程 （15）有，則（16） 假定懸浮泥沙的吸附量，對（11）式全斷面積分得： （17） 解決重金屬遷移轉(zhuǎn)換模型可將（12）至（17）式簡化成： 根據(jù)靜態(tài)吸附的結(jié)果，取b=0.543g/kg,k1=0.0076 1/(ppm·s)，k2=0.00084 1/s ,重金屬污染物的初始條件為：邊界條件為：在相同條件下，污染物的運(yùn)動遵循以下方程：。用向后顯式差分求解，結(jié)果如下圖2:由圖2可知：在不同時(shí)刻, 水相重金屬濃度由于泥沙

25、吸附, 峰值明顯減少, 而且由于泥沙吸附需要時(shí)間, 因此河段上泥沙吸附量的峰值總是滯后于水相濃度的峰值。從3中可以看出：污染物不受泥沙運(yùn)動影響, 在輸移過程中, 由于縱向綜合紊動擴(kuò)散作用, 峰型變化相同條件下重金屬污染物輸移, 盡管圖形相似, 但重金屬污染物的峰值因泥沙吸附而下降很多, 特別是后半部減小更為明顯, 且隨時(shí)間增加明顯加大。所以，重金屬在水體中的傳播特征是隨著離污染源的距離增加而減少。速度也隨之減小。2.土壤傳播模型的建立及求解： 重金屬在土壤中的傳播的影響因素也很多，不能單純的看作是線性變化，在此，我們基于重金屬在土壤中的質(zhì)量分?jǐn)?shù)，根據(jù)在水中傳播的方式類似的定義了一個(gè)非線性沉淀模

26、型。并推導(dǎo)出非線性沉淀情況下重金屬速度的控制方程。 對于飽和土壤空隙內(nèi)充填的溶液質(zhì)量分?jǐn)?shù)為，而土粒表面會吸附溶質(zhì)，因此，當(dāng)溶質(zhì)在土水界面達(dá)到平衡之后，會存在以下關(guān)系： （18）式中： 土壤中溶質(zhì)的總量在單位質(zhì)量土粒上的質(zhì)量分?jǐn)?shù)，包括土粒上吸附溶質(zhì)和孔隙溶液中的溶質(zhì)； 土顆粒表面所吸附的溶質(zhì)（j）的質(zhì)量分?jǐn)?shù)； 孔隙中的溶質(zhì)的質(zhì)量濃度（）； 土壤孔隙比； 土樣干密度（） 當(dāng)溶質(zhì)在土粒上的吸附呈線性狀況時(shí)，沉淀滿足Henry定律： （19）式中： 溶質(zhì)在土水界面的分配系數(shù)此時(shí)，滿足下式： （20）式中：、 均為模型參數(shù)，可表示為 土壤空隙率；S 土壤容重對于溶質(zhì)在土壤介質(zhì)上的吸附，其分配系數(shù)通常情況

27、下并非常數(shù)，而是隨溶質(zhì)濃度變化而變化的參數(shù)為了便于計(jì)算，本文對溶質(zhì)在土柱中土水間的分配關(guān)系做出如下簡化假定：， （21） 傳統(tǒng)的運(yùn)移控制方程如下： （22）式中：為時(shí)間，為阻滯因子，）表示隨土壤孔隙中溶質(zhì)的質(zhì)量濃度 而變化的非線性的擴(kuò)散系數(shù)，v為孔隙流體的流速式（22）是以土壤孔隙中溶質(zhì)的平衡濃度為基礎(chǔ)建立起來的假定擴(kuò)散系數(shù) 為常數(shù)，考慮吸附的非線性特性，將式（21）代入式（22），可得改進(jìn)的非線性重金屬運(yùn)移控制方程為 （23）其中： 基于土粒上溶質(zhì)質(zhì)量分?jǐn)?shù)的運(yùn)移參數(shù)提取方法，根據(jù)有限差分法，可將運(yùn)移方程（22）寫成如下差分格式： （23）式中： 方向被差分之后的差分網(wǎng)格節(jié)點(diǎn)的空間位置；，；、

28、為上下邊界處溶質(zhì)的質(zhì)量分?jǐn)?shù)；為0時(shí)刻土柱內(nèi)溶質(zhì)的質(zhì)量分?jǐn)?shù)上述差分格式可以利用Matlab軟件進(jìn)行求解，結(jié)合最小二乘法，利用內(nèi)嵌的最優(yōu)化工具箱求解運(yùn)移參數(shù)對比擬合結(jié)果與試驗(yàn)得到的濃度剖面數(shù)據(jù)，當(dāng)方差最小時(shí)即認(rèn)為此時(shí)的參數(shù)即為有代表性的典型參數(shù) 由模型可得重金屬在土壤中的傳播模型如下圖所示： 在圖中可以清晰地看出重金屬在土壤中的傳播隨時(shí)間的增減而減緩，而且速度較水體傳播慢。3. 大氣擴(kuò)散模型的建立及求解：大氣中重金屬主要來源于工業(yè)生產(chǎn)、汽車尾氣排放、汽車輪胎磨損產(chǎn)生的大量含有重金屬的有害氣體和粉塵。它們主要分布在工礦的周圍及公路、鐵路兩側(cè)。重金屬在大氣中的傳播一部分通過自然沉降和雨淋沉降進(jìn)入土壤圈

29、和水體中，一部分，直接污染空間。在這里因?yàn)閿?shù)據(jù)不足更為了簡化模型，我們忽略雨淋沉降，將通過大氣沉降地表的重金屬看做由自然沉降產(chǎn)生的。而自然沉降的主要作用是重力作用。在垂直方向，考慮粒子受到重力、浮力、流體阻力和隨機(jī)力的作用，建立粒子的運(yùn)動方程如下： （24）其中w為粒子垂直方向速度，為粒子受到的浮力，為粒子受到的流體阻力，為隨機(jī)力及大氣流動引起的隨機(jī)加速度，滿足如下條件: (25)粒子受到的阻力為，其中為動力粘性系數(shù)，r為粒子半徑。由式（24）得： 其中和分別為氣體和粒子密度，通常。令，由上式得到： 解此式可得： （26） 將 分解成為平均量和脈動量之和： 由這兩式可得： （27） （28）由

30、上式求出的方差，由于和不相關(guān)，則， （29） 由（28）式可得： 假設(shè)是平穩(wěn)隨機(jī)過程，則：由上面三式可得： 即：由（26）式可得： 講下式代入上式得： 將最后一項(xiàng)離散化得： （30） 令t=0得： 方差為： 與（29）式比較可得： 代入式（30）得到粒子運(yùn)動的離散方程： 在水平x、y方向不考慮重力作用，則在平均風(fēng)速為 (u，v，w)下粒子運(yùn)動的離散型方程為： 有題可知，結(jié)合計(jì)算精度和時(shí)間的考慮，時(shí)間步長 由公式： 。由于，在空氣中傳播的影響因素很多，而且其速度比水體快。結(jié)合上面三個(gè)模型可知，在同樣的條件下，金屬污染源在空氣中傳播最快，其次是在水中傳播，最后在土壤中傳播。而且這三種傳播狀態(tài)是相互

31、轉(zhuǎn)化的，對于同一污染源來說，對其進(jìn)行時(shí)間和空間的擬合得出如下傳播范圍（圖五）所示：由圖可知，污染源一定在三者疊加的正中心處，再結(jié)合八種重金屬元素在該區(qū)的分布圖，找出各金屬的污染源。如下圖（圖六）： 由以上八個(gè)圖可以看出，各暖色的區(qū)段為污染源分布區(qū)。其中這八種重金屬的污染源主要分布在工業(yè)區(qū)（坐標(biāo)為（3000,3000）附近）和交通區(qū)（坐標(biāo)為（13000,3000）附近），山區(qū)和公園綠地區(qū)幾乎沒有污染源。四、城市地質(zhì)環(huán)境的演變模式模型的建立：1.模型有缺點(diǎn)的分析：1）、優(yōu)點(diǎn)：a.在本模型中，運(yùn)用到了很多新思想。b.有許多問題先用理想化模型，之后再慢慢的增加約束條件，使問題更符合實(shí)際意義，更具有可推

32、廣性。c.在本文中，由于重金屬的傳播有固、液、氣三種相屬，而且，這三種相屬又是相互聯(lián)系，相互轉(zhuǎn)換的，不可很明確的區(qū)分開來，故，我們將三種模型非別建立，之后將其在空間和時(shí)間上擬合起來。使其更具說服力。2）、缺點(diǎn)：a.在本模型中，有很多假設(shè)，將模型理想化了，可能與實(shí)際有區(qū)別 b.個(gè)別數(shù)據(jù)題中未給出，我們經(jīng)過上網(wǎng)查詢，與實(shí)際問題會有誤差。2.城市地質(zhì)環(huán)境的演變模式模型： 經(jīng)過查詢可知：地質(zhì)環(huán)境的定義為：地殼上部包括巖石、水、氣和生物在內(nèi)的互相關(guān)聯(lián)的系統(tǒng)。本文前幾個(gè)模型主要研究城市表層土壤重金屬污染的問題，如果想要得到城市地質(zhì)環(huán)境的演變模式，僅有這些還遠(yuǎn)遠(yuǎn)不夠，因此需要收集一些重要信息才能更好地研究城

33、市地質(zhì)環(huán)境的演變模式。我們對地質(zhì)環(huán)境的主要影響因素進(jìn)行分析討論，得出如下（圖7）結(jié)果：地質(zhì)環(huán)境巖石圈水圈大氣圈生物圈地殼運(yùn) 動金屬含量生活工業(yè)污水地殼運(yùn)動廢氣污染紫外線植被凈化微生物凈化綠化程度 由圖可知，如果想研究城市地質(zhì)環(huán)境的演變模式，除了本題中的金屬影響之外，還有很多，如：地殼運(yùn)動、綠化程度、廢氣等等！因此，還需要收集一些地殼運(yùn)動、工業(yè)發(fā)展、綠化程度、紫外線輻射、城市人口的資料。同時(shí)還需要這個(gè)地區(qū)歷年的地質(zhì)環(huán)境影響因素資料及同一時(shí)期6個(gè)城區(qū)的相關(guān)資料。在擁有這些信息之后，我們通過以下方法建立并解決模型。 1.首先，用表示第i城區(qū)地質(zhì)環(huán)境的影響因素j的參數(shù)值。（i=1,2,6 j=1,2,

34、3,4）如下表所示：巖石圈參數(shù)水圈參數(shù)大氣圈參數(shù)生物圈參數(shù)1號城區(qū)y11y12y13y142號城區(qū)y21y22y23y243號城區(qū)y31y32y33y344號城區(qū)y41y42y43y445號城區(qū)y51y52y53y546號城區(qū)y61y62y63y64假設(shè)，四個(gè)參數(shù)所對應(yīng)的權(quán)重分別為，六個(gè)城區(qū)的地質(zhì)環(huán)境等級參數(shù)值記做。 在圖（7）中可以看出不能直接查出，也得經(jīng)過數(shù)據(jù)處理才能得到。但各影響因素的下屬因素能直接查出，同上方法可以通過已知信息計(jì)算出。 運(yùn)用公式 =* 可得各影響因素的權(quán)重值 2.之后，我們將這些權(quán)重值和一定時(shí)期影響因素參數(shù)值代入公式： （31） 進(jìn)行權(quán)重值的驗(yàn)算及調(diào)整。3.最后公式（3

35、1）和各已知數(shù)據(jù)算出城區(qū)不同時(shí)期地質(zhì)環(huán)境等級參數(shù)值。根據(jù)參數(shù)值就能夠研究出城市地質(zhì)環(huán)境的演變模式六模型的檢驗(yàn) 1. 本文中大多數(shù)圖形都用MATLAB 軟件繪制，而MATLAB軟件具有較強(qiáng)的仿真性，功能也十分強(qiáng)大，精確度也很高，從這個(gè)角度可證明我們結(jié)果的可靠性與方法的合理性。 2. 對于第一個(gè)重金屬的空間分布模型，將五個(gè)區(qū)域繪制在一幅圖表中，可以準(zhǔn)確地看出各個(gè)區(qū)域的聯(lián)系，有利于問題的分析。3. 模型一中計(jì)算重金屬的污染程度時(shí)，用平均數(shù)法和考慮重金屬毒性時(shí)的加權(quán)平均數(shù)法，二者對比得出相同結(jié)果，更說明了模型具有較強(qiáng)的說服力。而且，這兩種方法都比較簡單，便于理解。4. 在重金屬污染物傳播的模型中，本文

36、建立了三個(gè)傳播模型，而后將三種模型進(jìn)行時(shí)間和空間上的結(jié)合，定位污染源的位置，提高了模型的準(zhǔn)確性。5. 從模型與實(shí)際情況的對比來看，雖然我們建立模型有一些假設(shè)，但這些假設(shè)都是合乎常理的，重金屬的化學(xué)反應(yīng)是極其微弱的，故我們所建的模型的可靠性與真實(shí)性也是值得認(rèn)可的。七模型的推廣及發(fā)展前景文中第一個(gè)模型用到的平均數(shù)法和加權(quán)平均數(shù)法可以應(yīng)用于各個(gè)領(lǐng)域，如學(xué)校可以計(jì)算學(xué)生的學(xué)習(xí)成績的好壞，企業(yè)可以計(jì)算盈利和起伏的大小等。另外，本文中的重金屬傳播模型還可以推廣到消息的傳播、文化的傳播、疾病的傳播或者流言蜚語的傳播，從而對現(xiàn)今這個(gè)信息時(shí)代作一個(gè)正確的把握。尤其文中用MATLAB 軟件編程繪制的地形圖。用同樣

37、的方法可以運(yùn)用到地質(zhì)勘探和測量等方面。本文考慮的只是重金屬對土壤的污染問題，我們可以把它推廣到重金屬對植物和動物的影響，從而有利于對農(nóng)作物的培育和動物的養(yǎng)殖，甚至可以確定對人體帶來的危害，也可以應(yīng)用到其它金屬元素對土壤的污染和影響，從而研制促進(jìn)農(nóng)作物生長的化肥，有利于農(nóng)業(yè)的發(fā)展。此模型符合現(xiàn)今社會發(fā)展的狀況，目前，人類都關(guān)注著自己身體的健康，知道有害金屬對人體的危害，應(yīng)用本文的模型及分析方法，我們可以對某地區(qū)分析，選擇盡量遠(yuǎn)離易受重金屬污染的地理區(qū)域居住，盡量避免攝入含有重金屬元素的食物，有利于城市合理規(guī)劃，由此對人類的健康和安全帶來了可靠的理論和實(shí)踐依據(jù)。八模型的優(yōu)缺點(diǎn)一優(yōu)點(diǎn)： 1.本文中有

38、一些新的思想，如將三個(gè)模型在時(shí)間和空間上結(jié)合和了起來（也是合理解決此題的關(guān)鍵所在）。更準(zhǔn)確和客觀地找出了污染源的位置。2.文中采用的平均數(shù)法和加權(quán)平均數(shù)法比較簡單，使讀者可以正確體會其中的意思，并且，可以用初等數(shù)學(xué)的方法就可以解決。 3.對于重金屬的傳播問題，通過建立三個(gè)模型，進(jìn)行結(jié)合與比較，綜合考慮傳播的快慢和多少，準(zhǔn)確找出污染源的位置。4.研究城市地質(zhì)環(huán)境的模型中，考慮的因素比較多，更客觀全面地描繪了地質(zhì)環(huán)境的演變模式。二缺點(diǎn)： 1.本題中有些重金屬在水中還是有微弱的反應(yīng)的，故模型的建立和求解出來的答案與實(shí)際有一定的差距。2.由于題中沒有說明各個(gè)區(qū)域中的主要工商業(yè)，故在確定污染源的來源時(shí)不

39、能準(zhǔn)確定位。九參考文獻(xiàn)1 楊啟航,數(shù)學(xué)建模,北京:高等教育出版社,2003.2 韓中庚，數(shù)學(xué)建模方法及其應(yīng)用，北京：高等教育出版社，.2005.3 實(shí)用計(jì)算器4.50版軟件, 4 重金屬權(quán)數(shù)的確定，  dff04a1b0717fd5dd3a.5 解可新，最優(yōu)化方法， 天津：天津大學(xué)出版社，19966 實(shí)用MATLAB 軟件十、附件采樣點(diǎn)三維空間坐標(biāo)圖繪制程序d=data; % 只需從excel輸入三列數(shù)據(jù)，格式為： ，i表示行數(shù)x=d(:,1) % 采樣點(diǎn)坐標(biāo)x值y=d(:,2); % 采樣點(diǎn)坐

40、標(biāo)y值z=d(:,3); % 采樣點(diǎn)坐標(biāo)z(海拔)值nx=linspace(min(x),max(x),100);ny=linspace(min(y),max(y),100);xx,yy=meshgrid(nx,ny);zz=griddata(x,y,z,xx,yy,'v4');surfl(xx,yy,zz);shading interpcolormap(gray);hold on for i=1:319 for i=1:44 plot3(d(i,1),d(i,2),d(i,3),'ys'); end hold on for i=45:80 plot3(d(i,1),d(i,2),d(i,3),'y+'); end hold on for i=81:146 plot3(d(i,1),d(i,2),d(i,3),'bp'); end hold on for i=147:284 plot3(d(i,1),d(i,2),d(i,3),'ko'); end ho