獲獎?wù)撐?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)肅處理。我們參賽選擇的題號是（從A/B/C/D中選擇一項(xiàng)填寫）： A 我們的參賽報名號為（如果賽區(qū)設(shè)置

2、報名號的話）： 所屬學(xué)校（請?zhí)顚懲暾娜?太原工業(yè)學(xué)院 參賽隊(duì)員 (打印并簽名) ：1. 韓強(qiáng) 2. 鄭國軍 3. 曹玉乾 指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人 (打印并簽名)： 日期： 2011 年 9月 12日賽區(qū)評閱編號（由賽區(qū)組委會評閱前進(jìn)行編號）：2011高教社杯全國大學(xué)生數(shù)學(xué)建模競賽編 號 專 用 頁賽區(qū)評閱編號（由賽區(qū)組委會評閱前進(jìn)行編號）：賽區(qū)評閱記錄（可供賽區(qū)評閱時使用）：評閱人評分備注全國統(tǒng)一編號（由賽區(qū)組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進(jìn)行編號）：城市表層土壤重金屬污染分析摘要根據(jù)附件一中給出的海量數(shù)據(jù)，本文使用了matlab二維插值法建立八種重金屬三維

3、模型圖，通過模型圖得到8種重金屬元素在相同空間和相同元素在不同空間的分布。然后通過單因子污染指數(shù)法和內(nèi)梅羅綜合指數(shù)法，得出生活區(qū)和公園綠地區(qū)污染程度為輕污染；山區(qū)污染程度為尚清潔；工業(yè)區(qū)污染程度為重污染；主干道路區(qū)污染程度為中污染。本文使用excel軟件對數(shù)據(jù)進(jìn)行了處理與分析，得出城區(qū)污染的主要污染物是Hg和Pb。通過查閱大量的資料和文獻(xiàn)，分析得出重金屬污染的主要原因是制燒堿，汞化合物生產(chǎn)，做顏料，冶煉等工業(yè)“三廢”的排放。 本文分析了重金屬在表層土壤中的二維擴(kuò)散微分特征方程，由此建立了模型，確定出了污染源的位置，并且使用插值擬合的方法對模型進(jìn)行了檢驗(yàn)，確定出22號測量點(diǎn)所在的區(qū)域?yàn)橹饕廴驹?/p>

4、。本文所建立模型原理簡單，利用數(shù)據(jù)庫可方便地處理大量數(shù)據(jù)點(diǎn)；思路新穎，將每個離散點(diǎn)形象化的組織起來，適于研究區(qū)域內(nèi)離散數(shù)據(jù)。應(yīng)用廣泛，可在地學(xué)分析，計算機(jī)視覺，地理信息系統(tǒng)，計算幾何等領(lǐng)域進(jìn)行廣泛的研究和應(yīng)用。但在準(zhǔn)確定位（定位的具體方法）污染源或計算污染物的量以及確定二維偏微分方程的系數(shù)時考慮的因素較少，誤差略大。為此，本文使用因子分析，回歸分析，曲線擬合等方法建立了基本模型、考慮空氣模型和考慮空氣擴(kuò)散模型三種情形，給出了不同的解決方法，對原有模型進(jìn)行了改進(jìn)。關(guān)鍵詞：插值擬合法、二維擴(kuò)散微分方程模型、單因子污染指數(shù)法、回歸分析、質(zhì)量守恒原則 一、問題重述隨著城市經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷

5、增加，人類活動對城市環(huán)境質(zhì)量的影響日顯突出。對城市土壤地質(zhì)環(huán)境異常的查證，以及應(yīng)用查證獲得的海量數(shù)據(jù)資料開展城市環(huán)境質(zhì)量評價，研究人類活動影響下城市地質(zhì)環(huán)境的演變模式，日益成為人們關(guān)注的焦點(diǎn)。按照功能劃分，城區(qū)可分為生活區(qū)、工業(yè)區(qū)、山區(qū)、主干道路區(qū)及公園綠地區(qū)，將其記為1類區(qū)、2類區(qū)、5類區(qū)，不同的區(qū)域環(huán)境受人類活動影響的程度不同?，F(xiàn)對某城市城區(qū)土壤地質(zhì)環(huán)境進(jìn)行調(diào)查。為此，將所考察的城區(qū)劃分為間距1公里左右的網(wǎng)格子區(qū)域，按照每平方公里1個采樣點(diǎn)對表層土（010 厘米深度）進(jìn)行取樣、編號，并用GPS記錄采樣點(diǎn)的位置。應(yīng)用專門儀器測試分析，獲得了每個樣本所含的多種化學(xué)元素的濃度數(shù)據(jù)。另一方面，按照

6、2公里的間距在那些遠(yuǎn)離人群及工業(yè)活動的自然區(qū)取樣，將其作為該城區(qū)表層土壤中元素的背景值。附件1列出了采樣點(diǎn)的位置、海拔高度及其所屬功能區(qū)等信息，附件2列出了8種主要重金屬元素在采樣點(diǎn)處的濃度，附件3列出了8種主要重金屬元素的背景值?，F(xiàn)通過數(shù)學(xué)建模來完成以下任務(wù)：(1) 給出8種主要重金屬元素在該城區(qū)的空間分布，分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。(2) 通過數(shù)據(jù)分析，說明重金屬污染的主要原因。(3) 分析重金屬污染物的傳播特征，就此建立模型，確定污染源的位置。(4) 分析所建立模型的優(yōu)缺點(diǎn)，收集更多信息，更好地研究城市地質(zhì)環(huán)境的演變模式，建立跟完善的模型解決問題。二、問題分析這是一個關(guān)于分析

7、城區(qū)土壤重金屬污染分布和尋找污染源的問題。對于給出的海量數(shù)據(jù)，如何尋找數(shù)據(jù)間的關(guān)系和建立符合數(shù)據(jù)變化的函數(shù)模型是至關(guān)重要的，對此，我們使用各種數(shù)據(jù)處理軟件尋找數(shù)據(jù)間的關(guān)系，建立最優(yōu)的函數(shù)模型擬合數(shù)據(jù)的變化，進(jìn)而解決給出的問題。三、模型假設(shè)1. 所有附件數(shù)據(jù)真實(shí)可靠；2. 不考慮垂直方向的擴(kuò)散；3. 不考慮大氣、水流等其他因素對擴(kuò)散的影響；4. 不考慮植物對重金屬的吸收，表現(xiàn)為富集形式。四、定義與符號說明符號說明污染金屬指數(shù)污染金屬濃度實(shí)測值污染元素評價標(biāo)準(zhǔn)內(nèi)梅羅污染綜合指數(shù)評價因子個數(shù)元素污染指數(shù)的平均值元素污染指數(shù)的最大值i元素土壤污染累計指數(shù)i元素背景值五、模型的建立與求解第一部分：準(zhǔn)備工

8、作在進(jìn)行建模之前先對附件中的數(shù)據(jù)進(jìn)行了一些初步的處理。附件一數(shù)據(jù)的處理根據(jù)附件一中的海量數(shù)據(jù)，我們使用matlab軟件繪制出了五類區(qū)域和城區(qū)的地形圖（見附錄一）附件二數(shù)據(jù)的處理根據(jù)附件二中的數(shù)據(jù)，我們使用matlab軟件，繪制出了各元素在各區(qū)的分布情況圖，以及各元素在全城區(qū)的分布情況圖。a各元素在各區(qū)的情況分布圖，附圖如下：（見附錄二）b各元素在城區(qū)情況分布圖：（見附錄三） 第二部分：問題一的解答（一）問題1模型的建立和求解1.為了得出8種重金屬元素在該城區(qū)的空間分布，我們基于對附件一及附件二中海量數(shù)據(jù)的分析，以及查證到的信息，利用matlab軟件，繪制出了8種重金屬元素在該城區(qū)的不同區(qū)域的分

9、布情況圖（見附錄一），以及某元素在全城區(qū)的分布情況圖（見附錄二）。a、由8種重金屬元素在該城區(qū)的不同區(qū)域的分布情況圖，我們可以清晰的了解到各區(qū)的重金屬分布情況一類區(qū)中，相對于其他重金屬元素，Zn的分布量最多，其次為Cd，Cr，Pb，Hg，Cu，Ni，As。在二類區(qū)中，相對于其他重金屬元素，Hg的分布量明顯高于其他重金屬元素，Cu的分布量居其次，之后為Zn，Cd，Pb，Cr，Ni，As。三類區(qū)中，相對于其他重金屬元素，Pb的分布量略高于其他重金屬元素，其他元素分布量基本相同。在四類區(qū)中，相對于其他重金屬元素，Hg的分布量出奇的高，其次為Zn，Cd，Cu，cr，Pb，Ni，As。在第五類中，相對于

10、其他重金屬元素，Hg，Zn，Cd的含量高于其他元素，Pb，Cu，Cr遠(yuǎn)低于前三種元素，Ni，As最少。b、由單個元素在全城區(qū)的分布情況圖，我們可以清晰的了解到全城區(qū)各種金屬元素的分布情況在城區(qū)中，就污染量而言，Hg高居榜首，Zn居其后，之后為Cu，Cd，Cr，Pb，Ni，As。就污染范圍而言，Hg，Zn，Cu，Cd的范圍很大，Cr，Ni，Pb，As略小。2.我們利用現(xiàn)有資料以及查找到的信息，利用單因子污染指數(shù)法和內(nèi)梅羅綜合指數(shù)法來評價城區(qū)內(nèi)不同區(qū)域金屬的污染程度。a單因子污染指數(shù)法是國內(nèi)通用的一種重金屬污染評價方法，其計算公式如下：我們采用土壤環(huán)境質(zhì)量標(biāo)準(zhǔn)（GB15618-1995）中國家一級

11、標(biāo)準(zhǔn)作為環(huán)境質(zhì)量的標(biāo)準(zhǔn)（見表1-1-1）。 含量單位：項(xiàng)目一級二級三級土壤PH值自然背景<6.56.57.5>7.5>7.5As 1530252030Cd 0.20.30.30.61Cr 90250300350400Cu 3550100100400Hg 0.150.30.511.5Ni 40405060200Pb 35250300350500Zn 100200250300500 表1-1-1<1表明未受污染，>1表示已經(jīng)受到污染，數(shù)值越大，表明受到的污染越嚴(yán)重。利用附錄二中的數(shù)據(jù)，我們使用Excel軟件做出了城區(qū)內(nèi)不同區(qū)域內(nèi)重金屬污染物的單因子污染指數(shù)（見表1-1

12、-2）。單因子污染指數(shù)和內(nèi)梅羅污染綜合指數(shù)表單因子污染指數(shù)一類區(qū)二類區(qū)三類區(qū)四類區(qū)五類區(qū)As 0.41803030.4834260.2696060.38 0.42 Cd 1.449806821.9655560.76 1.81 1.402714Cr 0.766871210.59 0.43 0.6461580.484844Cu 1.41 3.64 0.4947791.7862570.86262Hg 0.62 4.2823690.273042.9876160.766611Ni 0.458556820.4952920.3863450.4403210.382243Pb1.974467532.658311

13、.0444551.8211351.734531Zn2.370086362.7792750.7329422.4409941.54最大單因子污染指數(shù)2.370086364.2823691.0444552.9876161.734531平均單因子污染指數(shù)1.183701222.1126930.5494561.5386970.949196內(nèi)梅羅污染綜合指數(shù)1.873293623.3765480.8345022.3762831.40 表1-1-2根據(jù)表1-1-2的數(shù)據(jù)，我們利用Exce軟件l做出8種重金屬元素在五類區(qū)中的單因子污染指數(shù)柱狀圖，如下圖： 由上圖可見，在一類區(qū)中，Cd、Cu、Pb、Zn的單因子

14、污染指數(shù)均高于1,表明一類區(qū)已受污染；在二類區(qū)中，Cd、Pb、Cu、Hg、Zn的單因子污染指數(shù)均高于1，表明二類區(qū)已經(jīng)受到污染；在三類區(qū)中，僅Pb的單因子污染指數(shù)略高于1，其為受污染區(qū)；在四類區(qū)中，Cd、Pb、Cu、Hg、Zn的單因子污染指數(shù)均高于1，表明其為受污染區(qū)；在五類區(qū)中，Cd、Pb、Zn的單因子污染指數(shù)均高于1，所以五類區(qū)均為受污染區(qū)。綜合以上分析，可以得出，城區(qū)的不同區(qū)域均受到不同程度的污染。b單因子指數(shù)可以判斷出環(huán)境的主要污染因子，但環(huán)境是一個復(fù)雜的體系，環(huán)境污染往往由多個污染因子復(fù)合污染導(dǎo)致，故單因子指數(shù)法不能全面、綜合的反應(yīng)土壤的污染程度，基于此，我們采用內(nèi)梅羅污染綜合指數(shù)法

15、進(jìn)行分析。內(nèi)梅羅綜合指數(shù)法計算公式如下：根據(jù)附件中的數(shù)據(jù)及單因子污染指數(shù)，我們用Excel做出了城區(qū)不同區(qū)域的重金屬污染內(nèi)梅羅綜合指數(shù)，見表1-1-2據(jù)表1-1-2我們用Excel繪制出了重金屬內(nèi)梅羅污染綜合指數(shù)，如下圖： 圖1土壤污染指數(shù)分級內(nèi)梅羅綜合污染指數(shù)全面反映了各污染物對土壤污染的不同程度，同時充分考慮高含量物質(zhì)對本土壤的影響，結(jié)合相關(guān)標(biāo)準(zhǔn)，將土壤污染等級進(jìn)行劃分（見表1-1-3）。等級P污染程度污染水平10.7安全清潔20.71.0警戒線尚清潔31.02.0輕污染超過背景值，視輕污染42.03.0中污染土壤開始受中度污染53.0重污染污染已相當(dāng)嚴(yán)重 表1-1-3從上圖1及表1-1-

16、3，可以得出三類區(qū)的內(nèi)梅羅污染指數(shù)剛過警戒線，尚清潔；一類區(qū)和五類區(qū)的內(nèi)梅羅污染指數(shù)介于1.02.0，屬于輕污染；四類區(qū)的內(nèi)梅羅污染指數(shù)介于2.03.0，屬中污染；二類區(qū)的內(nèi)梅羅污染指數(shù)高于3.0，屬于重污染。c由單因子污染指數(shù)法和內(nèi)梅羅綜合指數(shù)法得出的結(jié)論進(jìn)行綜合的分析可以得出：二類區(qū)，四類區(qū)污染程度比其他三區(qū)嚴(yán)重，尤以二類區(qū)最為嚴(yán)重；一類區(qū)，五類區(qū)屬輕污染；三類區(qū)污染程度最輕。第三部分：問題二模型的建立與求解1根據(jù)附錄中所給數(shù)據(jù)及查證獲得的數(shù)據(jù)信息，我們首先計算出了土壤污染累積指數(shù)，見表2-1-1：土壤累積污染指數(shù)計算公式：全區(qū)平均值背景平均值標(biāo)準(zhǔn)偏差范圍污染物累積指數(shù)As (g/g)5.

17、903.60.91.85.41.64Cd (ng/g)295.4913030701902.27Cr (g/g)52.6431913491.70Cu (g/g)57.3913.23.66.020.44.35Hg (ng/g)267.9035819517.65Ni (g/g)17.3012.33.84.719.91.41Pb (g/g)64.6331619432.08Zn (g/g)197.31691441972.86表2-2-1各區(qū)重金屬比值圖 根據(jù)上表中給出的污染累積指數(shù)，及上圖我們可以推斷，城區(qū)的重金屬污染，與人類的生產(chǎn)、生活有密切關(guān)系，通過上圖將全區(qū)金屬濃度平均值與背景值做比較，可以初步得

18、出此城區(qū)的重金屬污染物主要為Cd、Cu、Hg、Pb、Zn?；诖耍摮菂^(qū)的污染可能為工業(yè)排放和交通尾氣排放。2.基于對問題一中所得數(shù)據(jù)的分析，我們繪制出了各元素在各區(qū)的單因子污染指數(shù)柱狀圖，見下圖：據(jù)上圖分析，該城區(qū)不同功能區(qū)土壤中重金屬污染貢獻(xiàn)不同，從圖中可以清晰的得出，Pb的污染范圍最廣，基本對各個區(qū)都有影響，Cd、Zn影響4個區(qū)， Cu影響3個區(qū)，Hg影響2個區(qū)，As、Cr、Ni對各個區(qū)基本沒有影響。 由問題一得出的結(jié)論繪制圖表如下：各區(qū)單因子污染指數(shù)分布圖各區(qū)主要污染元素污染貢獻(xiàn)元素主要污染貢獻(xiàn)元素一類區(qū)Cd、Cu、Pb、ZnZn二類區(qū)Cd、Pb、Cu、Hg、ZnHg三類區(qū)PbPb四類

19、區(qū)Cd、Pb、Cu、Hg、ZnHg五類區(qū)Cd、Pb、ZnPb據(jù)上圖和表分析，該城區(qū)不同功能區(qū)土壤中重金屬污染貢獻(xiàn)具有一定的差別，在一類區(qū)中，Zn、Pb的污染貢獻(xiàn)最大，Cd、Cu居其次，根據(jù)查閱到的文獻(xiàn)資料，Cd、Cu、Pb、Zn來源主要為工業(yè)廢水的排放和農(nóng)藥的使用，由于一類區(qū)屬于生活區(qū)，故可判斷其污染的主要原因可能為生產(chǎn)活動中化肥施用過量，或是城市垃圾的焚燒，燃煤等。在二類區(qū)中，Hg、Cu、Zn、Pb的污染貢獻(xiàn)最大，Cd居次，Hg、Cu、Zn、Pb的主要來源于工業(yè)廢水，污泥，農(nóng)藥等，由于二類區(qū)為工業(yè)區(qū)，可基本判斷為工業(yè)污染，主要原因可能為顏料廠，冶金廠，電鍍廠等工業(yè)廢水、廢渣處理不合格，亂排亂

20、放。在三類區(qū)中，僅Pb具有一定的污染貢獻(xiàn)，且其含量很低，與背景值最為接近，鑒于三類區(qū)為山區(qū)，其Pb污染可能為氣流運(yùn)動，沉降累計。在四類區(qū)中，Cd、Pb、Cu、Hg、Zn的污染貢獻(xiàn)普遍很大，由于四類區(qū)為主干道路區(qū)，其污染原因最為可能是汽車尾氣的排放，以及汽車輪胎的磨損。在五類區(qū)中，Cd、Pb、Zn的污染貢獻(xiàn)大，五類區(qū)為公園綠地區(qū)，其污染原因最為可能是通過氣流運(yùn)動，降水等沉降積累。綜合以上分析，該城區(qū)的重金屬污染主要原因在于工業(yè)釋放及交通排放。第四部分：問題三的重金屬傳播特征的模型1模型的建立和求解此模型的主要目的是根據(jù)測得的污染數(shù)據(jù),得到一個關(guān)于重金屬污染物傳播的偏微分方程,進(jìn)而確定污染源位置的

21、計算方法。我們根據(jù)給出的各個觀測點(diǎn)的位置坐標(biāo)和污染數(shù)據(jù)利用對重金屬污染的線性插值，畫出了重金屬污染圖，將附件三給出的的八種重金屬元素的值作為重金屬污染物的背景濃度（），觀測點(diǎn)的處理如下：第一，光滑插值后如果有不光滑數(shù)據(jù)出現(xiàn)，考慮到計算機(jī)處理數(shù)據(jù)的誤差，予以忽略。第二，分別對八種元素處理，得到各元素的污染坐標(biāo)，對于相對誤差不超過5%的坐標(biāo)予以合并，取其平均值。第三，單個元素分析時，只考慮權(quán)重達(dá)到0.8以上的污染數(shù)據(jù)，其他數(shù)據(jù)因其影響太小，予以忽略?？紤]到污染物在土壤中的運(yùn)動狀態(tài)，我們忽略了重金屬污染的垂直擴(kuò)散，假設(shè)重金屬均勻擴(kuò)散且是一維的，以污染源為原點(diǎn)，觀測坐標(biāo)X軸方向作為X軸，設(shè)污染物在原點(diǎn)

22、的污染為連續(xù)過程，選用一維場中污染物質(zhì)的二維擴(kuò)散微分方程為：初始條件為：邊界條件為：當(dāng)處瞬時點(diǎn)源的解析解：其中是hantush函數(shù)，且常量參數(shù)：以污染源位置、 、為變量，結(jié)合格點(diǎn)搜索算法，反復(fù)調(diào)整變量的值，直到收斂函數(shù)的值最小，其中為個測量點(diǎn)的污染值；為背景值。2.模型的求解此模型預(yù)測結(jié)果表明：As的污染源在第84個和第178個測量點(diǎn)附近；Cd的污染源在第22個和第95個測量點(diǎn)附近；Cr的污染源在第22個測量點(diǎn)附近；Cu的污染源在第8個和第22個測量點(diǎn)附近；Hg的污染源在第9個和第182個和第257個測量點(diǎn)附近；Ni的污染源在第22個測量點(diǎn)附近；Pb的污染源在第6個和第16個測量點(diǎn)附近；Zn的

23、污染源在第61個測量點(diǎn)附近；分析結(jié)果：可看出：主要污染來自于工業(yè)區(qū)和交通區(qū)，且在交通區(qū)測得的第22個測量點(diǎn)污染最嚴(yán)重，說明22號附近的區(qū)域是主要污染源，與前文分析結(jié)果一致，說明模型正確。3.模型驗(yàn)證：我們使用matlab軟件的griddata命令下的'cubic'（以三角形為基礎(chǔ)的三次方程內(nèi)插）插值的方式對海量數(shù)據(jù)進(jìn)行了處理和繪圖，得到了各元素的分布圖，經(jīng)過matlab的計算，可得到污染最嚴(yán)重處的坐標(biāo)，從圖中的位置可判斷污染源的位置如下：As分布圖：可得和兩個點(diǎn)，對比前面的模型，誤差為Cd分布圖：可得和兩個點(diǎn)，對比前面的模型，誤差為Cr分布圖：可得一個點(diǎn)，對比前面的模型，誤差為

24、Cu分布圖：可得和兩個點(diǎn)，對比前面的模型，誤差為Hg分布圖：可得 、和三個點(diǎn)，對比前面的模型，誤差為Ni分布圖：可得一個點(diǎn)，對比前面的模型，誤差為Pb分布圖：可得和兩個點(diǎn)，對比前面的模型，誤差為Zn分布圖：可得一個點(diǎn)，對比前面的模型，誤差為由此可看出，模型的準(zhǔn)確性稍差。第五部分：問題四的求解模型的優(yōu)缺點(diǎn)模型的優(yōu)點(diǎn)：模型原理簡單，利用數(shù)據(jù)庫可方便地整理大量數(shù)據(jù)點(diǎn)。模型思路新穎，將每個離散點(diǎn)形象化的組織起來，適于研究區(qū)域內(nèi)離散數(shù)據(jù)。應(yīng)用廣泛，可在地學(xué)分析，計算機(jī)視覺，地理信息系統(tǒng)，計算幾何等領(lǐng)域進(jìn)行廣泛的研究和應(yīng)用。模型的缺點(diǎn)：不能準(zhǔn)確定位（定位的具體方法）污染源或計算污染物的量。確定二維偏微分方

25、程的系數(shù)時考慮的因素較少，誤差略大。加入的信息鑒于以上缺點(diǎn)，我們認(rèn)為，還應(yīng)收集以下信息來完善模型的準(zhǔn)確性，1.重金屬在土壤擴(kuò)散過程中的揮發(fā)因素。2.如果在重金屬污染源附近有河流流經(jīng)，應(yīng)考慮河流的沖刷作用。3.大氣的各種現(xiàn)象（刮風(fēng)，下雨，沙塵暴等）也會影響重金屬的擴(kuò)散。4.整個城區(qū)的土地中是否有富集的因素影響。5.土壤的粘滯力也會對重金屬的擴(kuò)散產(chǎn)生重要影響。模型的修改：根據(jù)質(zhì)量守恒原則，以污染物在土壤中的擴(kuò)散為主要影響因素建立了具有恒定擴(kuò)散系數(shù)的二維擴(kuò)散方程，該模型為三題模型的改進(jìn)，仍然假設(shè)擴(kuò)散僅發(fā)生在土壤的表面。式中，C為土壤中某污染物的濃度；Z為土壤中測定點(diǎn)到土壤表面的距離；D為污染物在土壤

26、中的擴(kuò)散系數(shù)；t為時間。對于不同的邊界條件，上式有不同的解。a邊界條件當(dāng)t=0且時，C=；當(dāng)Z=0且t>0時，C=0；當(dāng)Z=L時，此時，擴(kuò)散方程的解為：揮發(fā)量為：式中，L為土壤層的厚度；其余符號同前。上式中，當(dāng)指數(shù)項(xiàng)>10時，該式的累加項(xiàng)將遠(yuǎn)小于1，可以忽略，則上式可以簡化為：這時，t時間內(nèi)單位面積上由污染物的揮發(fā)總量為：對于土壤中重金屬污染物的揮發(fā)，某時刻、某位置處的濃度為：式中，erf為誤差函數(shù)符號，上式適用的范圍為b邊界條件當(dāng)t=0且時，C=；當(dāng)t=0且Z>0時，C=0；當(dāng)t=0且Z=0時，C=0此時，土壤中的重金屬污染物的濃度計算模型如下：化學(xué)污染物的揮發(fā)量為:當(dāng)?shù)闹?/p>

27、很大時，上式同樣可化簡為的形式。c如果假設(shè)重金屬污染物在土壤深度L之內(nèi)是均勻變化的，且在下部邊界沒有污染物的流失，而土壤中污染物的揮發(fā)速率受土壤上部空氣去處揮發(fā)出得污染物速率的控制，即相對于流失而言，擴(kuò)散的影響可以忽略；如果空氣停止流動，則污染物的揮發(fā)量為零?；谏鲜黾僭O(shè)，有以下邊界條件：當(dāng)t=0且時，C=；當(dāng)Z=L時，；當(dāng)t>0且Z=0時，偏微方程在此邊界條件下的解為：式中，v為土壤上放空氣流動速度，；為土壤上放空氣中重金屬污染物的濃度，；為吸附等溫系數(shù)，為重金屬污染物在空氣中于土壤中的濃度比，即；為土壤中濃度；為下式的根多數(shù)情況下，有，可簡化為：上式中，為下式的根：重金屬污染物通過土

28、壤表面的揮發(fā)量為：如果土壤上方空氣流動不明顯，揮發(fā)出的重金屬污染物在土壤上方穩(wěn)定空氣層中的擴(kuò)散系數(shù)不能忽略，對于這種情況，可通過設(shè)定下述邊界條件球的其計算模式。邊界條件：當(dāng)t=0且時，C=；當(dāng)Z=L時，；當(dāng)t>0且Z=0時，C=；當(dāng)t>0且Z=0時， 解為：式中，為下式的根：重金屬污染物通過土壤表面的揮發(fā)量為： 上述模型是建立在質(zhì)量平衡原則基礎(chǔ)上的，是從宏觀的角度描述污染物從土壤中的揮發(fā)規(guī)律的。六、模型的評價與推廣本模型采用了多種方法有效的減少計算量，使計算的時間減少到可以接受的程度。模型原理簡單可方便的整理大量數(shù)據(jù)點(diǎn)，能良好的處理區(qū)域內(nèi)離散數(shù)據(jù)問題，在地學(xué)分析、計算機(jī)視覺、地理信

29、息系統(tǒng)、計算機(jī)和等方面有著廣泛的應(yīng)用。但是模型仍然可以進(jìn)一步優(yōu)化，由于時間倉促，水平有限，沒能將這些想法完全實(shí)現(xiàn)。七、參考文獻(xiàn)1 鄭彤,陳春云，環(huán)境系統(tǒng)數(shù)學(xué)模型，北京：化學(xué)工業(yè)出版社，2003。2 殷劍宏，吳開亞，圖論及其算法，合肥：中國科學(xué)技術(shù)大學(xué)出版社，2004。3 朱旭，李煥琴，籍萬新，MATLAB軟件與基礎(chǔ)數(shù)學(xué)實(shí)驗(yàn)，西安：西安交通大學(xué)出版社，2008。4 方影，孫慶文，高等數(shù)學(xué)與數(shù)學(xué)模型，北京：高等教育出版社，2009。5鄔學(xué)軍，周凱，宋軍全，數(shù)學(xué)建模競賽輔導(dǎo)教程，杭州：浙江大學(xué)出版社，2009。6 趙睿新，環(huán)境污染化學(xué)，北京：化學(xué)工業(yè)出版社，2004。7 胡克林，張芬榮，呂貽忠，王茹

30、，徐艷，北京市大興區(qū)土壤重金屬含量的空間分布特征，環(huán)境科學(xué)學(xué)報，第24卷第3期：464468,2005年4月。 8段雪梅，蔡煥興，巢文軍，南京市表層土壤重金屬污染特征及污染來源，環(huán)境科學(xué)與管理，第35卷第17期：3177,2010年10月。9國家環(huán)境保護(hù)局，土壤環(huán)境質(zhì)量標(biāo)準(zhǔn)，北京：中國環(huán)境出版社，1995。10情系文檔心，09年全國數(shù)學(xué)建模優(yōu)秀論文，八、附錄附錄一 城區(qū)地形示意圖如下：附錄二 一區(qū)各元素的分布圖：1-1-1 1-1-2 二區(qū)各元素分布圖： 2-1-12-1-2 三區(qū)各元素分布圖：3-1-13-1-2 四區(qū)各元素分布圖：4-1-14-1-2 五區(qū)各元素分布圖：5-1-15-1-2

31、附錄三 As在城區(qū)情況分布圖：1-2-1 Cd在城區(qū)情況分布圖：2-2-1 Cr在城區(qū)情況分布圖：3-2-1 Cu在城區(qū)情況分布圖：4-2-1 Hg在城區(qū)情況分布圖：5-2-1 Ni在城區(qū)情況分布圖：6-2-1 Pb在城區(qū)情況分布圖：7-2-1 Zn在城區(qū)情況分布圖：8-2-1附錄四%1區(qū)As分布圖x=data(1:44,2);y=data(1:44,3);z=data(1:44,6);ti=0:10:26453;tu=0:10:18449;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi

32、)hold onplot3(x,y,z,'o')%1區(qū)Cd分布圖x=data(1:44,2);y=data(1:44,3);z=data(1:44,7);ti=0:10:26453;tu=0:10:18449;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%1區(qū)Cr分布圖x=data(1:44,2);y=data(1:44,3);z=data(1:44,8);ti=0:10:26453;tu=0:10:184

33、49;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%1區(qū)Cu分布圖x=data(1:44,2);y=data(1:44,3);z=data(1:44,9);ti=0:10:26453;tu=0:10:18449;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o'

34、;)%1區(qū)Hg分布圖x=data(1:44,2);y=data(1:44,3);z=data(1:44,10);ti=0:10:26453;tu=0:10:18449;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%1區(qū)Ni分布圖x=data(1:44,2);y=data(1:44,3);z=data(1:44,11);ti=0:10:26453;tu=0:10:18449;xi,yi=meshgrid(ti,tu); zi

35、=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%1區(qū)Pb分布圖x=data(1:44,2);y=data(1:44,3);z=data(1:44,12);ti=0:10:26453;tu=0:10:18449;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%1區(qū)Zn分布圖x=data(1:44,2);y=

36、data(1:44,3);z=data(1:44,13);ti=0:10:26453;tu=0:10:18449;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%2區(qū)As分布圖x=data(45:80,2);y=data(45:80,3);z=data(45:80,6);ti=0:10:22674;tu=0:10:17980;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,

37、9;cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%2區(qū)Cd分布圖x=data(45:80,2);y=data(45:80,3);z=data(45:80,7);ti=0:10:22674;tu=0:10:17980;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%2區(qū)Cr分布圖x=data(45:80,2);y=data(45:80,3);z=data(4

38、5:80,8);ti=0:10:22674;tu=0:10:17980;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%2區(qū)Cu分布圖x=data(45:80,2);y=data(45:80,3);z=data(45:80,9);ti=0:10:22674;tu=0:10:17980;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi

39、,yi,zi)hold onplot3(x,y,z,'o')%2區(qū)Hg分布圖x=data(45:80,2);y=data(45:80,3);z=data(45:80,10);ti=0:10:22674;tu=0:10:17980;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%2區(qū)Ni分布圖x=data(45:80,2);y=data(45:80,3);z=data(45:80,11);ti=0:10:226

40、74;tu=0:10:17980;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%2區(qū)Pb分布圖x=data(45:80,2);y=data(45:80,3);z=data(45:80,12);ti=0:10:22674;tu=0:10:17980;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot

41、3(x,y,z,'o')%2區(qū)Zn分布圖x=data(45:80,2);y=data(45:80,3);z=data(45:80,13);ti=0:10:22674;tu=0:10:17980;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%3區(qū)As分布圖x=data(81:146,2);y=data(81:146,3);z=data(81:146,6);ti=0:10:28654;tu=0:10:18397

42、;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%3區(qū)Cd分布圖x=data(81:146,2);y=data(81:146,3);z=data(81:146,7);ti=0:10:28654;tu=0:10:18397;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o

43、')%3區(qū)Cr分布圖x=data(81:146,2);y=data(81:146,3);z=data(81:146,8);ti=0:10:28654;tu=0:10:18397;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%3區(qū)Cu分布圖x=data(81:146,2);y=data(81:146,3);z=data(81:146,9);ti=0:10:28654;tu=0:10:18397;xi,yi=meshg

44、rid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%3區(qū)Hg分布圖x=data(81:146,2);y=data(81:146,3);z=data(81:146,10);ti=0:10:28654;tu=0:10:18397;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%3區(qū)Ni

45、分布圖x=data(81:146,2);y=data(81:146,3);z=data(81:146,11);ti=0:10:28654;tu=0:10:18397;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%3區(qū)Pb分布圖x=data(81:146,2);y=data(81:146,3);z=data(81:146,12);ti=0:10:28654;tu=0:10:18397;xi,yi=meshgrid(ti,tu

46、); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%3區(qū)Zn分布圖x=data(81:146,2);y=data(81:146,3);z=data(81:146,13);ti=0:10:28654;tu=0:10:18397;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%4區(qū)As分布圖x=data

47、(147:284,2);y=data(147:284,3);z=data(147:284,6);ti=0:10:26424;tu=0:10:17949;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%4區(qū)Cd分布圖x=data(147:284,2);y=data(147:284,3);z=data(147:284,7);ti=0:10:26424;tu=0:10:17949;xi,yi=meshgrid(ti,tu); zi

48、=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%4區(qū)Cr分布圖x=data(147:284,2);y=data(147:284,3);z=data(147:284,8);ti=0:10:26424;tu=0:10:17949;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%4區(qū)Cu分布圖x=data(14

49、7:284,2);y=data(147:284,3);z=data(147:284,9);ti=0:10:26424;tu=0:10:17949;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%4區(qū)Hg分布圖x=data(147:284,2);y=data(147:284,3);z=data(147:284,10);ti=0:10:26424;tu=0:10:17949;xi,yi=meshgrid(ti,tu); zi=g

50、riddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%4區(qū)Ni分布圖x=data(147:284,2);y=data(147:284,3);z=data(147:284,11);ti=0:10:26424;tu=0:10:17949;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%4區(qū)Pb分布圖x=data(147

51、:284,2);y=data(147:284,3);z=data(147:284,12);ti=0:10:26424;tu=0:10:17949;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%4區(qū)Zn分布圖x=data(147:284,2);y=data(147:284,3);z=data(147:284,13);ti=0:10:26424;tu=0:10:17949;xi,yi=meshgrid(ti,tu); zi=griddata(x,y,z,xi,yi,'cubic');mesh(xi,yi,zi)hold onplot3(x,y,z,'o')%5區(qū)As分布圖x=data(285:319,2);y=data(285:319,3);z=data(285:319,6);ti=0:10:2