城市表層土壤重金屬污染分析數(shù)學(xué)建模國(guó)家二等獎(jiǎng)

2011高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽承諾書(shū)我們仔細(xì)閱讀了中國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽的競(jìng)賽規(guī)則.我們完全明白，在競(jìng)賽開(kāi)始后參賽隊(duì)員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊(duì)外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問(wèn)題。我們知道，抄襲別人的成果是違反競(jìng)賽規(guī)則的,如果引用別人的成果或其他公開(kāi)的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻(xiàn)的表述方式在正文引用處和參考文獻(xiàn)中明確列出。我們鄭重承諾，嚴(yán)格遵守競(jìng)賽規(guī)則，以保證競(jìng)賽的公正、公平性。如有違反競(jìng)賽規(guī)則的行為，我們將受到嚴(yán)肅處理。我們參賽選擇的題號(hào)是（從A/B/C/D中選擇一項(xiàng)填寫(xiě)）： 我們的參賽報(bào)名號(hào)為（如果賽區(qū)設(shè)置報(bào)名號(hào)的話）：所屬學(xué)校（請(qǐng)?zhí)顚?xiě)完整的全名）：參賽隊(duì)員(打印并簽名)：1.2.3.指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人(打印并簽名)：日期：年月日賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：2011高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽編號(hào)專用頁(yè)賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：賽區(qū)評(píng)閱記錄（可供賽區(qū)評(píng)閱時(shí)使用）：評(píng)閱人評(píng)分備注全國(guó)統(tǒng)一編號(hào)（由賽區(qū)組委會(huì)送交全國(guó)前編號(hào)）：全國(guó)評(píng)閱編號(hào)（由全國(guó)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：城市表層土壤重金屬污染分析摘要本文以某城區(qū)為例，對(duì)其土壤地質(zhì)環(huán)境重金屬含量進(jìn)行取樣分析，找出取樣點(diǎn)被污染的主要原因并尋找到取樣點(diǎn)的污染源。為此，對(duì)于第（1）和第（2）小題我們對(duì)取樣點(diǎn)的數(shù)據(jù)統(tǒng)計(jì)分析，用MATLAB軟件畫(huà)出8張?jiān)氐姆植紙D，各代表同一重金屬元素在不同功能區(qū)不同含量范圍內(nèi)的分布點(diǎn)，統(tǒng)計(jì)出點(diǎn)的數(shù)目，對(duì)均值和方差進(jìn)行對(duì)比得出該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度和主要污染原因。針對(duì)第（3）小題，我們建立了2個(gè)模型來(lái)研究，干塵擴(kuò)散模型——模型一，高斯水流擴(kuò)散模型——模型二。模型一是為了求解重金屬元素附著在干塵粒上的擴(kuò)散方式和范圍，運(yùn)用煙霧擴(kuò)散模型，結(jié)論是根據(jù)取樣點(diǎn)的某金屬含量程度來(lái)判斷該點(diǎn)距污染源的距離，則污染源在以該距離為半徑的圓的圓周上，在根據(jù)取樣點(diǎn)的金屬含量高低分布情況，找到污染源。模型二是研究重金屬在土壤中的擴(kuò)散情況，在相間物質(zhì)交換為平衡的條件下,可用阻滯系數(shù)來(lái)表示其影響，可得到一個(gè)微分方程帶入?yún)?shù)，可以找到某金屬的污染源坐標(biāo)。關(guān)鍵詞：煙霧擴(kuò)散模型圖例統(tǒng)計(jì)污染物運(yùn)移物質(zhì)交換MATLAB1問(wèn)題的提出隨著城市經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷增加，人類活動(dòng)對(duì)城市環(huán)境質(zhì)量的影響日顯突出。對(duì)城市土壤地質(zhì)環(huán)境異常的查證，以及如何應(yīng)用查證獲得的海量數(shù)據(jù)資料開(kāi)展城市環(huán)境質(zhì)量評(píng)價(jià)，研究人類活動(dòng)影響下城市地質(zhì)環(huán)境的演變模式，日益成為人們關(guān)注的焦點(diǎn)。按照功能劃分，城區(qū)一般可分為生活區(qū)、工業(yè)區(qū)、山區(qū)、主干道路區(qū)及公園綠地區(qū)等，分別記為1類區(qū)、2類區(qū)、……、5類區(qū)，不同的區(qū)域環(huán)境受人類活動(dòng)影響的程度不同。本文根據(jù)某個(gè)城市的5個(gè)區(qū)域的采樣數(shù)據(jù)分析出8種主要重金屬元素在該城區(qū)的空間分布，并分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。通過(guò)數(shù)據(jù)分析，說(shuō)明重金屬污染的主要原因。分析重金屬污染物的傳播特征，由此建立模型，確定污染源的位置。附件1列出了采樣點(diǎn)的位置、海拔高度及其所屬功能區(qū)等信息，附件2列出了8種主要重金屬元素在采樣點(diǎn)處的濃度，附件3列出了8種主要重金屬元素的背景值。現(xiàn)要求通過(guò)數(shù)學(xué)建模來(lái)完成以下任務(wù)：(1)給出8種主要重金屬元素在該城區(qū)的空間分布，并分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。(2)通過(guò)數(shù)據(jù)分析，說(shuō)明重金屬污染的主要原因。(3)分析重金屬污染物的傳播特征，由此建立模型，確定污染源的位置。(4)分析你所建立模型的優(yōu)缺點(diǎn)，為更好地研究城市地質(zhì)環(huán)境的演變模式，還應(yīng)收集什么信息？有了這些信息，如何建立模型解決問(wèn)題？2問(wèn)題的分析對(duì)于問(wèn)題和問(wèn)題，主要采用圖形對(duì)照法，根據(jù)附件中的數(shù)據(jù)使用MATLAB圖形處理軟件圖形化，分析出同一重金屬元素在不同功能區(qū)及不同金屬元素在同一功能區(qū)的含量分布情況，找出污染的主要原因；而對(duì)于問(wèn)題，為了更好的研究重金屬傳播特性和確定污染源位置，我們建立了兩種模型：模擬重金屬在空氣中的傳播，建立干塵擴(kuò)散模型——模型一；考慮到金屬在土壤里的傳播，建立土壤擴(kuò)散模型——模型二。在問(wèn)題中，對(duì)模型的優(yōu)缺點(diǎn)進(jìn)行了評(píng)價(jià)，提出了其它在模型中沒(méi)有考慮到的一些會(huì)影響重金屬污染程度的因素，并在模型一和模型二基礎(chǔ)上進(jìn)一步優(yōu)化。問(wèn)題一：金屬元素符號(hào)說(shuō)明：砷鎘鉻銅汞鎳鉛鋅重金屬含量分段說(shuō)明：、、、、、、、含量安全范圍分別為：、、、、、、、，我們分段時(shí)從均值開(kāi)始到安全范圍上限，以后的范圍斷是人為分段（由于相關(guān)數(shù)據(jù)難以查到），具體分布可以參照表1經(jīng)過(guò)參考附錄中的坐標(biāo)點(diǎn)及海拔我們用軟件實(shí)現(xiàn)了該城區(qū)的三維模擬再根據(jù)已知數(shù)據(jù)用畫(huà)出俯瞰的功能區(qū)分布圖如圖2示，經(jīng)過(guò)圖形處理，得到該市大概的地區(qū)分布圖形如示：由砷的相關(guān)數(shù)據(jù)畫(huà)出各個(gè)區(qū)砷含量分布：由和知砷在部分工業(yè)區(qū)含量在以上，屬于重度污染，還有部分工業(yè)區(qū)砷含量在，交通區(qū)砷含量普遍在內(nèi)，屬于輕度污染；由鎘的相關(guān)數(shù)據(jù)得出該城市鎘的含量分布用比，可以得出該城市在山區(qū)鎘含量為，屬于輕度污染；在公園和生活區(qū)污染相對(duì)要嚴(yán)重，含鎘量在；在交通區(qū)含鎘量在；而工業(yè)區(qū)含鎘量則在之間，污染相當(dāng)嚴(yán)重。由鉻元素的數(shù)據(jù)可以得到鉻元素的分布：通過(guò)和對(duì)比，該城市普遍鉻元素含量在，在生活區(qū)含量在，而工業(yè)區(qū)和交通區(qū)污染則更嚴(yán)重，含量在之間；由附件中銅元素的數(shù)據(jù)可以畫(huà)出銅元素在該城市的分布：通過(guò)和的對(duì)比，發(fā)現(xiàn)在工業(yè)區(qū)和交通區(qū)銅元素含量在之間，其它功能區(qū)含量普遍低于；根據(jù)附件中汞元素的分布數(shù)據(jù)得出：與對(duì)比后可以看出工業(yè)區(qū)﹑交通區(qū)以及地處工業(yè)區(qū)附近的公園綠地區(qū)含量偏高，在之間，生活區(qū)和山區(qū)汞含量則低于；從附件中鎳元素的數(shù)據(jù)得到了鎳元素在該城區(qū)的分布：通過(guò)和比較，工業(yè)區(qū)和交通區(qū)以及生活區(qū)鎳含量偏高，在之間，而其它功能區(qū)則普遍低于；同樣的，由附件中的數(shù)據(jù)可以得到鉛和鋅元素在城市的分布和：通過(guò)和圖3對(duì)比鉛以及鋅在交通區(qū)﹑業(yè)區(qū)以及生活區(qū)含量高，鉛含量在80μg/g～480μg/g之間，鋅含量則是500μg/g～3800μg/g的范圍內(nèi)；上述圖對(duì)比為宏觀上的判斷，針對(duì)此，我們列出了同一重金屬元素在不同含量段在整個(gè)城區(qū)的分布，并統(tǒng)計(jì)出每個(gè)含量段的取樣點(diǎn)數(shù)，如表1As(μg/g)Cd(ng/g)Cr(μg/g)Cu(μg/g)含量范圍樣點(diǎn)數(shù)含量范圍樣點(diǎn)數(shù)含量范圍樣點(diǎn)數(shù)含量范圍樣點(diǎn)數(shù)3.6-5.488130-1906431-4914813.2-20.4585.4-1513190-50015949-1009920.4-15021115-30.135500-100039100-2005150-500101000-20007200-5003500-10000500-100021000-26002Hg(ng/g)Ni(μg/g)Pb(μg/g)Zn(μg/g)含量范圍樣點(diǎn)數(shù)含量范圍樣點(diǎn)數(shù)含量范圍樣點(diǎn)數(shù)含量范圍樣點(diǎn)數(shù)35-516112.3-19.916831-438669-977151-1509719.9-507043-8010997-16068150-5004150-100180-15053160-50094500-10007100-1501150-2509500-1500131000-16003250-48031500-38005表1通過(guò)表1的分析，、、、、在該城區(qū)屬于重度污染；、達(dá)到污染；分析，、、、、在該城區(qū)屬于重度污染；、達(dá)到污染；屬于輕度污染綜上，得出重金屬分布情況：八種重金屬元素在工業(yè)區(qū)、交通區(qū)含量較高，污染較嚴(yán)重；在生活區(qū)和公園綠地區(qū)含量相對(duì)要少，屬于輕度污染；而山區(qū)各種金屬元素含量除了接近工業(yè)區(qū)和交通區(qū)的區(qū)域，重金屬含量都在安全范圍內(nèi)。問(wèn)題2：主要思路：分別對(duì)五個(gè)功能區(qū)重金屬元素的分布情況進(jìn)行分析，找出該區(qū)某重金屬含量超過(guò)安全范圍的主要原因。通過(guò)均值和偏差的對(duì)比得討論如下：功能區(qū)均值砷（As）鎘（Cd）鉻（Cr）銅（Cu）汞（Hg）鎳（Ni）鉛（Pb）鋅（Zn）生活區(qū)7.48502.97118.16174.8953.6623.37116.22501.39工業(yè)區(qū)4.21323.8350.7950.4778.9814.6960.96264.86山區(qū)5.11224.9539.5925.0890.5214.9444.81104.51交通區(qū)5.83286.7243.9641.61312.8117.7157.92162.75公園綠地區(qū)5.39236.0538.9118.4147.4614.8741.0592.31標(biāo)準(zhǔn)偏差砷（As）鎘（Cd）鉻（Cr）銅（Cu）汞（Hg）鎳（Ni）鉛（Pb）鋅（Zn）生活區(qū)3.93281.54166.97415.523081.3519.4996.14547.46工業(yè)區(qū)1.48162.9745.2452.79146.886.1225.11607.79山區(qū)3.71249.7422.1321.75261.599.7933.31165.03交通區(qū)2.46173.1614.7829.891719.685.3531.21159.96公園綠地區(qū)2.39189.9016.338.0448.645.8812.4267.64對(duì)照?qǐng)D和表所列數(shù)據(jù)，得以下的結(jié)論：1、生活區(qū)鉛和鋅的含量相對(duì)其他重金屬含量較高，考慮到人們的生活環(huán)境，如人們所用到的家具和房間的涂料，以及家用汽車的增加導(dǎo)致尾氣的排放、生活垃圾等都含有大量的鉛；含鋅化合物主要用作油漆、油墨、塑料、廣告顏料、化妝品、蓄電池的原料或填充劑。所以可以看出該市的人們?cè)谏顓^(qū)使用的涂料、化妝品、電池和相關(guān)商品較多并且處理不的不完善，導(dǎo)致含鋅量較高。2、工業(yè)區(qū)所有重金屬的含量都很高。長(zhǎng)期以來(lái)我國(guó)由于科技和經(jīng)濟(jì)欠發(fā)達(dá)，居民的環(huán)境保護(hù)意識(shí)薄弱，工業(yè)生產(chǎn)過(guò)程中產(chǎn)生的廢水、廢氣、廢渣未經(jīng)處理直接排放。其中含有大量的重金屬，經(jīng)過(guò)自然的沉降雨水的淋溶等途徑進(jìn)入土壤導(dǎo)致Hg、Pb、Cu、Cr、As、Ni、Zn、Cd含量逐年增加導(dǎo)致工業(yè)區(qū)所有重金屬含量都很高。3、山區(qū)各重金屬含量大部分都在正常范圍內(nèi)，部分地區(qū)屬于輕度污染。主要原因是受到工業(yè)區(qū)和生活區(qū)的影響。工業(yè)區(qū)通過(guò)空氣傳播和水流傳播以及生活區(qū)的生活垃圾污染導(dǎo)致部分山區(qū)的輕度污染。4、交通區(qū)鉛和鎘、銅、鋅的含量較高。交通區(qū)汽車燃油、潤(rùn)滑油和鍍金部分的燃燒或磨損釋放出大量重金屬元素。一般以道路為中心成帶狀分布，導(dǎo)致路邊土壤中這幾種重金屬元素明顯高于非交通線路兩側(cè)的土壤。5、公園綠地區(qū)汞、鋅、鉛的含量相對(duì)較高。首先考慮到該區(qū)所處的地理位置,地處生活區(qū)、工業(yè)區(qū)、交通區(qū)附近。所以公園較嚴(yán)重的重金屬污染可能來(lái)源于交通、工業(yè)污染和居民生活廢棄物等污染。問(wèn)題3：一模型的假設(shè)模型一干塵擴(kuò)散模型1、假設(shè)無(wú)風(fēng)情況下，重金屬看作在地面某一點(diǎn)向四周等強(qiáng)度的在無(wú)窮空間擴(kuò)散，不計(jì)地形影響。2、重金屬物質(zhì)附著在干沉上的擴(kuò)散服從擴(kuò)散定律，即單位時(shí)間通過(guò)單位法平面面積的流量與它的濃度梯度成正比。3、重金屬物質(zhì)在穿過(guò)降雨時(shí)，其強(qiáng)度由于雨水的吸收而減少，減少比率為常數(shù)。4、風(fēng)的影響只考慮重金屬污染源的某一風(fēng)向、某一風(fēng)速的影響。風(fēng)的影響范圍，隨著時(shí)間而減小。模型二土壤水流擴(kuò)散模型1、重金屬在土壤多孔介質(zhì)中有三種存在形式，假設(shè)相間物質(zhì)交換為平衡2、延滯系數(shù)唯一。二、模型的建立1.模型概述重金屬擴(kuò)散的范圍和強(qiáng)度主要受到污染源的自然擴(kuò)散、風(fēng)以及降雨的影響。自然擴(kuò)散符合煙霧擴(kuò)散模型，在給定污染源總含量等數(shù)據(jù)時(shí)，可以求出靜態(tài)擴(kuò)散范圍。風(fēng)力具有加大擴(kuò)散面積的作用，而降雨具有吸收稀釋重金屬元素含量的作用。2.參數(shù)示意污染源濃度污染源重金屬擴(kuò)散的速度表示某點(diǎn)的濃度污染源重金屬擴(kuò)散的總量擴(kuò)散的時(shí)間降雨對(duì)重金屬的吸收率有風(fēng)時(shí)候的速度有風(fēng)時(shí)候的方向風(fēng)力影響隨時(shí)間減小量泄漏點(diǎn)距有效地面的高度(m)分別為用濃度標(biāo)準(zhǔn)差表示的軸上的擴(kuò)散參數(shù)3.模型的求解求解主要思想：在無(wú)降雨情況下污染源重金屬總量等于擴(kuò)散后的總量，所以有，當(dāng)有降雨的情況下，擴(kuò)散后金屬含量有變化，吸收后剩余量為.3.1模型一建立過(guò)程建立空間坐標(biāo)系：將重金屬開(kāi)始擴(kuò)散時(shí)設(shè)為，以污染源為坐標(biāo)原點(diǎn)，建立以正東為x正方向，正北為y正方向，建立三維直角坐標(biāo)系，如下：3.1.1推導(dǎo)重金屬隨干沉擴(kuò)散范圍半徑（無(wú)風(fēng)無(wú)降雨的擴(kuò)散模型）時(shí)刻t無(wú)窮空間中任意一點(diǎn)的重金屬含量記為P。根據(jù)假設(shè)2，單位時(shí)間通過(guò)單位法相面的流量為：·····················s是擴(kuò)散系數(shù)，grad表示梯度，負(fù)號(hào)表示由濃度高向濃度低的地方擴(kuò)散。考察空間域Ω，Ω的體積為V，包圍Ω的曲面為S，S的外法線向量為n,則在[t,t+Δt]內(nèi)通過(guò)Ω的流量為：··················而Ω內(nèi)重金屬物質(zhì)含量的增量為：·····················由擴(kuò)散前后金屬總量不變得：································根據(jù)曲面面積分的奧氏公式：·······················其中是散度記號(hào)。由可得：無(wú)界區(qū)域的拋物型偏微分方程······················根據(jù)假設(shè)1，初始條件為作用在坐標(biāo)原點(diǎn)的點(diǎn)源函數(shù)，可以記作：···························得：·······這個(gè)結(jié)果表明。對(duì)于任意時(shí)刻t，重金屬物質(zhì)濃度C的等值面是球面x2+y2+z2=R2,并且隨著球面半徑R的增加C的只是連續(xù)減少的；當(dāng)R—>∞或t—>∞時(shí)，C—>0通過(guò)將x2+y2+z2=R2代入方程（8）可以解出R：··············R為在給定擴(kuò)散系數(shù)s，影響時(shí)間t，重金屬含量為C，以及污染源總含量W的情況下的受自然擴(kuò)散影響范圍的半徑，即在不考慮風(fēng)力、降雨情況下的影響范圍。通過(guò)matlab，可以實(shí)現(xiàn)該模型的圖像生成，如圖一。在此假定：s=0.00001(擴(kuò)散系數(shù))，W=3000（μg/g）（污染源總量，假定），C=103.20（重金屬在檢測(cè)點(diǎn)物質(zhì)濃度，這里取鉛），t=1000000(秒)（擴(kuò)散時(shí)間）r=0.5(假設(shè)雨水對(duì)于重金屬吸收率為50%)ws=0.04(公里每小時(shí))（假設(shè)作用在擴(kuò)散源的風(fēng)速為40米/小時(shí)）wa=0.785（假設(shè)風(fēng)向與x0y面平行，且與x軸成正45度角）wd=0.3（假設(shè)風(fēng)力影響范圍從R開(kāi)始，每天減弱300米）我們調(diào)用函數(shù)ZJ_Shu(k,Q,C,t,r,wv,we,wj)（函數(shù)本身可參考附錄），加入?yún)?shù)為ZJ_Shu(0.0001,3000,103.2,1000000,0.5,0.04,0.785,0.3)，符號(hào)意思看參數(shù)說(shuō)明表，結(jié)果為：經(jīng)過(guò)去不同重金屬的均值，我們得到元素取值（平均值）半徑RAs3.624.9564×10^3Cd13015.3267×10^3Cr3118.1306×10^3Cu13.220.4128×10^3Hg3517.8477×10^3Ni12.320.6261×10^3Pb3118.1306×10^3Zn6916.4310×10^33.1.2考慮風(fēng)和降雨影響的模型建立降雨因素：假設(shè)降雨對(duì)于放射性物質(zhì)的吸收率為r,則方程（9）變?yōu)椋猴L(fēng)力因素：只考慮產(chǎn)生在輻射原點(diǎn)某一角度，某一強(qiáng)度的風(fēng)力影響。風(fēng)力的影響范圍從R開(kāi)始每天減少，并且對(duì)于影響范圍內(nèi)每一點(diǎn)的作用相同。3.2模型結(jié)果分析模型的分析結(jié)果說(shuō)明在給定的條件下，擴(kuò)散的最遠(yuǎn)距離為24.95公里。由于現(xiàn)實(shí)中風(fēng)向的不確定性，出于安全考慮，在污染源周邊24.95公里范圍的地域內(nèi)都有污染，而根據(jù)距離，可以判斷出污染源在以距離為半徑的圓周上，尋找到工業(yè)區(qū)或者是生活區(qū)。模型二：重金屬元素在土壤中的傳播模型分析：重金屬元素在土壤多孔介質(zhì)中有三種可能的存在形式：:溶于水中、揮發(fā)為氣體及吸附于固體顆粒。揮發(fā)性有機(jī)污染物在水、氣、固體顆粒三相間的物質(zhì)交換與分配是決定其運(yùn)移的重要因素,在相間物質(zhì)交換為平衡的條件下,可用阻滯系數(shù)來(lái)表示其影響。污染物與土壤固體顆粒間的非平衡吸附解吸是相間交換中影響污染物運(yùn)移的最重要的機(jī)制,由于天然土壤具有固有的不均質(zhì)性,必須用多個(gè)反應(yīng)系數(shù)才能準(zhǔn)確描述污染物與固體顆粒間的非平衡吸附解吸。在土壤環(huán)境中,一系列的機(jī)制控制著污染物的運(yùn)移:一、地下水流決定了污染物的運(yùn)動(dòng)方向和速率;二、擴(kuò)散使污染物產(chǎn)生縱向及橫向的轉(zhuǎn)移;三、污染物與土壤顆粒中有機(jī)質(zhì)及礦物質(zhì)之間的吸附解吸、污染物在土壤包氣帶中的水氣界面處的物質(zhì)交換使污染物的運(yùn)移受到阻滯作用;四、由于具有揮發(fā)性,污染物還隨氣體遷移、擴(kuò)散;五、土壤中的生物與化學(xué)作用使污染物降解,生成無(wú)害物質(zhì)或其他有害物質(zhì)。要預(yù)測(cè)污染物的運(yùn)移和其歸宿,必須對(duì)土壤—污染物—空氣—水這一復(fù)雜的相互作用的系統(tǒng)及污染物在其中的諸多遷移機(jī)制有充分理解,并把該系統(tǒng)模型化從而得到污染源的位置。模型的建立符號(hào)說(shuō)明氣體中的污染物濃度亨利系數(shù)水中的污染物濃度顆粒表面的污染物濃度污染物在水與固體顆粒間的分配系數(shù)分配系數(shù)土壤含水率土壤含氣率土壤干容重孔隙水流速度水平距離擴(kuò)散系數(shù)比例系數(shù)距離傳播速度時(shí)間某已知點(diǎn)的橫坐標(biāo)某已知點(diǎn)的縱坐標(biāo)污染源的橫坐標(biāo)污染源的縱坐標(biāo)1氣—水—固三相間的物質(zhì)交換在有水氣共存的土壤中,即土壤非飽和帶中,污染物隨水運(yùn)移時(shí),水中的污染物會(huì)與氣體及土壤顆粒發(fā)生物質(zhì)交換。水氣兩相間的物質(zhì)交換平衡可用亨利定律來(lái)表示式中,為氣體中的污染物濃度,為水中的污染物濃度,H為亨利系數(shù)。不同物質(zhì)的H值可在文獻(xiàn)中查到。對(duì)于水中污染物與固體顆粒間的物質(zhì)交換,經(jīng)常用到以下這一簡(jiǎn)單的線性方程:=式中,為固體顆粒表面的污染物濃度,為污染物在水與固體顆粒間的分配系數(shù)。有機(jī)物質(zhì)在水與固體顆粒間的分配與土壤顆粒中的有機(jī)質(zhì)含量、土壤顆粒的大小與表面積等因素直接相關(guān),而在天然土壤中,土壤的這些性質(zhì)往往不是均一的,用線性的方程來(lái)描述水—固間的物質(zhì)分配與實(shí)測(cè)數(shù)據(jù)往往不很吻合,而非線性方程則更好地表示了其間的物質(zhì)交換。方程是一個(gè)經(jīng)常被使用的非線性吸附等溫線方程,其表達(dá)式為:式中,為分配系數(shù),表示水—固吸附非線性指數(shù)。方程、都屬經(jīng)驗(yàn)公式,其中的參數(shù)值可通過(guò)專門(mén)的實(shí)驗(yàn)測(cè)出。污染物由氣相吸附在固體顆粒表面的量一般很小,可以忽略不計(jì)。假設(shè)氣相為靜止,且不考慮污染物在氣相中的擴(kuò)散,污染物在非飽和帶中隨水運(yùn)移可由下列一維對(duì)流—擴(kuò)散方程來(lái)表示:式中,，分別為土壤含水率與含氣率,為土壤干容重,為孔隙水流速度,為水平距離,為擴(kuò)散系數(shù)。方程左邊第二、第三項(xiàng)分別表示污染物在固相與氣相中的分配。在水—?dú)?、水—固之間的物質(zhì)分配為平衡的條件下,可把方程、(或)代入,且方程兩邊同除以得:為延滯系數(shù),延滯系數(shù)Rt代表污染物運(yùn)移時(shí)由于發(fā)生了水—?dú)饧八腆w顆粒之間的物質(zhì)交換而產(chǎn)生的延滯作用。在不同的土壤含水率條件及不同的污染物運(yùn)移承載相(氣或水)情況下,延滯系數(shù)的表達(dá)式是不同的,如考慮污染物隨氣相在非飽和帶中運(yùn)移,則可由污染物通過(guò)土柱的時(shí)間與不發(fā)生水—固或氣—固吸附時(shí)非揮發(fā)性試劑通過(guò)土柱的時(shí)間的對(duì)比實(shí)驗(yàn)來(lái)求得。近來(lái)有實(shí)驗(yàn)表明,實(shí)驗(yàn)測(cè)得的污染物在非飽和土柱中隨氣相運(yùn)移時(shí)的延滯系數(shù)大于由求得的預(yù)測(cè)值。這部分額外的延滯效應(yīng)則應(yīng)歸于污染物在水—?dú)饨缑嫔系姆e累。污染物在水—?dú)饨缑嫔系姆e累不僅使污染物運(yùn)移受阻,而且其影響還相當(dāng)大,如等在實(shí)驗(yàn)中發(fā)現(xiàn)(三氯乙烯)的滯后是由于在水—?dú)饨缑嫣幍姆e累造成的。因此,揮發(fā)性有機(jī)污染物在水—?dú)狻倘嗟姆峙淇赡馨ㄒ韵聨讉€(gè)方面:污染物在水—?dú)鈨上嚅g的分配、污染物在水—固體顆粒兩相間的分配及污染物在水—?dú)饨缑嫣幍睦鄯e。同時(shí)考慮這三個(gè)方面后的污染物隨氣體運(yùn)移的延滯系數(shù)為:式中,為污染物在水氣界面處累積的吸附系數(shù),為有效水氣界面面積。由所求出的延滯系數(shù)可以看出和即濃度有關(guān)，我們可以假設(shè)延滯系數(shù)即擴(kuò)展系數(shù)，顯然延滯系數(shù)和傳播速度成正比而從以污染源為起點(diǎn)延滯系數(shù)將逐漸降低，設(shè)污染源的位置為任意一地點(diǎn)（采樣的點(diǎn)）位置為。單位時(shí)間傳播的距離的積分就是污染源到這點(diǎn)的距離從而由得根據(jù)和已知數(shù)據(jù)便可以求出污染源的位置。用MATLAB求解并取代表性坐標(biāo)如元素坐標(biāo)As(μg/g)（5013，7202）Cd(ng/g)（2724，2576）（4034，6005）Cr(μg/g)（4632，4315）（3504，6013）Cu(μg/g)（2367，3443）（3845，6012）Hg(ng/g)（4812，1845）（4865，7011）Ni(μg/g)（4055，6089）Pb(μg/g)（2002，2951）（4990，5310）Zn(μg/g)（14024，9529）（3970，5978）（5372，7537）經(jīng)過(guò)數(shù)據(jù)分析，主要污染源在該城區(qū)的西南方向，集中了工業(yè)區(qū)，交通區(qū)，生活區(qū)，海拔在4——8之間。問(wèn)題44模型的評(píng)價(jià)對(duì)于干塵擴(kuò)散模型——模型一優(yōu)點(diǎn):考慮了無(wú)風(fēng)無(wú)雨時(shí)候的理想情況和有風(fēng)有雨的情況，把現(xiàn)實(shí)生活中的大部分情況都包括在內(nèi)，并且得出了不同情況時(shí)候的傳播關(guān)系。缺點(diǎn)：考慮的因素較少，忽略了人為，濕度，自然環(huán)境的吸收能力等因素。得到的結(jié)果的準(zhǔn)確性自然較低。對(duì)于土壤擴(kuò)散模型——模型二優(yōu)點(diǎn)：通過(guò)文獻(xiàn)，得到了直觀且簡(jiǎn)單的計(jì)算方法和關(guān)系式。缺點(diǎn)：也只簡(jiǎn)單的求出延滯系數(shù)，通過(guò)所求的式子得到各個(gè)量之間的關(guān)系，得到濃度和位置的一個(gè)關(guān)系式。對(duì)不同的重金屬的傳播途徑的分析和利用不夠透徹。為更好地研究城市地質(zhì)環(huán)境的演變模式，還應(yīng)收集的信息有:本地區(qū)不同地方的降水量、地基的復(fù)雜程度、人們的活動(dòng)范圍以及將來(lái)社會(huì)的發(fā)展趨向（如哪里新建住房等）等。有了以上這些信息。我們認(rèn)為應(yīng)綜合的根據(jù)以下各方面的原因分析。（1）降雨——降雨的多少必定影響這個(gè)地區(qū)的土壤和污染等；（2）地基的復(fù)雜程度——地基的不同也對(duì)這個(gè)地區(qū)的承受能力和變化有影響；（3）人們的活動(dòng)范圍——人們的活動(dòng)對(duì)這個(gè)地方的演變起到重要的作用，人多的地方演變的必定更快。（4）將來(lái)社會(huì)的發(fā)展趨向——社會(huì)的發(fā)展趨向?qū)Νh(huán)境的影響更是至關(guān)重要，直接影響人們活動(dòng)的地域，從而影響環(huán)境。5模型的推廣1、本模型可以擴(kuò)展到其他污染物的污染，只需檢測(cè)出相應(yīng)相應(yīng)的污染物在各個(gè)區(qū)域的含量分布，就可以找出污染源；2、本模型對(duì)擴(kuò)散型污染物都有一定的適用性，可以參照此解決一些擴(kuò)散性污染的問(wèn)題；3、本模型具有較強(qiáng)的應(yīng)用價(jià)值。環(huán)境的污染容易破壞生態(tài)鏈，從而影響人類的生活，重金屬的污染，甚至可以直接危害人類的健康。該模型可以通過(guò)控制一些因素找出被污染城市的污染源，從而及時(shí)制定應(yīng)急措施，減小損失。參考文獻(xiàn)楊桂元，數(shù)學(xué)模型應(yīng)用實(shí)例，合肥：工業(yè)大學(xué)出版社，2007年6月第一版汪國(guó)強(qiáng)，數(shù)學(xué)建模優(yōu)秀案例選編，華南：華南理工大學(xué)出版社，2001年8月第一版周博謝東來(lái)，MATLAB科學(xué)計(jì)算，三河市：機(jī)械工業(yè)出版社，2010年5月第一版附錄：圖形實(shí)現(xiàn)主代碼A=[......];%此矩陣數(shù)據(jù)太多，在此省略bianhao=A(:,1);%編號(hào)As=[A(:,2)];Cd=[A(:,3)];Cr=[A(:,4)];Cu=[A(:,5)];Hg=[A(:,6)];Ni=[A(:,7)];Pb=[A(:,8)];Zn=[A(:,9)];x=[A(:,10)];y=[A(:,11)];z=[A(:,12)];gnqu=[A(:,13)];X_hang=length(gnqu);%scatter(x,y,5,z)%散點(diǎn)圖[X,Y,Z]=griddata(x,y,z,linspace(0,28654)',linspace(0,18449),'v4');%插值%pcolor(X,Y,Z);%shadinginterp%偽彩色圖%所有監(jiān)測(cè)點(diǎn)網(wǎng)絡(luò)圖fig1=figure(1);set(fig1,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高線圖holdonplot(x,y,'-*')title('檢測(cè)點(diǎn)網(wǎng)絡(luò)圖');xlabel('x');ylabel('y');fig2=figure(2);set(fig2,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高線圖holdon%Ask=0;fori=1:X_hangif3.6<=As(i)&As(i)<=5.4plot(x(i),y(i)+250,'^');holdonelseif5.4<As(i)&A(i)<=15plot(x(i)+250,y(i),'s','markerfacecolor','g');holdonelseifAs(i)>15plot(x(i),y(i)-250,'d','markerfacecolor','r');text(x(i)+200,y(i)+200,sprintf('(%0.1f)',As(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendtitle('As的分布','backgroundcolor','y','color','b','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');text(1000,17900,'

-重度污染As(μg/g)>15','backgroundcolor','k','color','w');text(1000,16900,'△-輕度污染3.6<=As(μg/g)<=5.4','backgroundcolor','k','color','w');text(1000,15900,'□-污染5.4<As(μg/g)<=15','backgroundcolor','k','color','w');%Cdfig3=figure(3);set(fig3,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高線圖holdonfori=1:X_hangifCd(i)<130&Cd(i)<190plot(x(i),y(i),'v','markerfacecolor','k');%text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cd(i)),'color','k','fontsize',10,'fontweight','bold');holdonelseif190<=Cd(i)&Cd(i)<=500plot(x(i),y(i),'p','markerfacecolor','g','markersize',10);%text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cd(i)),'color','k','fontsize',10,'fontweight','bold');holdonelseif500<Cd(i)&Cd(i)<=1000plot(x(i),y(i),'o','markerfacecolor','r');%text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cd(i)),'color','r','fontsize',10,'fontweight','bold');holdonelseif1000<Cd(i)&Cd(i)<2000plot(x(i),y(i),'d','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cd(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendtitle('Cd的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');text(100,17900,'▲-130<Cd(ng/g)<190','color','w','backgroundcolor','k');text(100,16900,'★-190<=Cd(ng/g)<=500','color','w','backgroundcolor','k');text(100,15900,'●-500<Cd(ng/g)<=1000','color','w','backgroundcolor','k');text(100,14900,'◆-1000<Cd(ng/g)<2000','color','w','backgroundcolor','k');%Crfig4=figure(4);set(fig4,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高線圖holdonfori=1:X_hangif31<=Cr(i)&Cr(i)<=49plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif49<Cr(i)&Cr(i)<=100plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif100<Cr(i)&Cr(i)<=200plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif200<Cr(i)&Cr(i)<=500plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif500<Cr(i)&Cr(i)<=1000plot(x(i),y(i),'s','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cr(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendendtext(100,17900,'▲-32<Cr(μg/g)<49','color','w','backgroundcolor','k');text(100,16900,'★-49<=Cr(μg/g)<=100','color','w','backgroundcolor','k');text(100,15900,'●-100<Cr(μg/g)<=200','color','w','backgroundcolor','k');text(100,14900,'◆-200<Cr(μg/g)<500','color','w','backgroundcolor','k');text(100,13900,'■-500<Cr(μg/g)<1000','color','w','backgroundcolor','k');title('Cr的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Cufig5=figure(5);set(fig5,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高線圖holdonfori=1:X_hangif13.2<=Cu(i)&Cu(i)<=20.4plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif20.4<Cu(i)&Cu(i)<=150plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif150<Cu(i)&Cu(i)<=500plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif500<Cu(i)&Cu(i)<=1000plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif1000<Cu(i)&Cu(i)<=2600plot(x(i),y(i),'s','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cu(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendendtext(100,17900,'▲-13.2<Cu(μg/g)<20.4','color','w','backgroundcolor','k');text(100,16900,'★-20.4<=Cu(μg/g)<=150','color','w','backgroundcolor','k');text(100,15900,'●-150<Cu(μg/g)<=500','color','w','backgroundcolor','k');text(100,14900,'◆-500<Cu(μg/g)<1000','color','w','backgroundcolor','k');text(100,13900,'■-1000<Cu(μg/g)<2600','color','w','backgroundcolor','k');title('Cu的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Hgfig6=figure(6);set(fig6,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高線圖holdonfori=1:X_hangif35<=Hg(i)&Hg(i)<=51plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif51<Hg(i)&Hg(i)<=150plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif150<Hg(i)&Hg(i)<=500plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif500<Hg(i)&Hg(i)<=1000plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif1000<Hg(i)&Hg(i)<=1600plot(x(i),y(i),'s','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Hg(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendendtext(100,17900,'▲-35<Hg(ng/g)<51','color','w','backgroundcolor','k');text(100,16900,'★-51<=Hg(ng/g)<=150','color','w','backgroundcolor','k');text(100,15900,'●-150<Hg(ng/g)<=500','color','w','backgroundcolor','k');text(100,14900,'◆-500<Hg(ng/g)<1000','color','w','backgroundcolor','k');text(100,13900,'■-1000<Hg(ng/g)<=1600','color','w','backgroundcolor','k');title('Hg的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Nifig7=figure(7);set(fig7,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高線圖holdonfori=1:X_hangif12.3<=Ni(i)&Ni(i)<=19.9plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif19.9<Ni(i)&Ni(i)<=50plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif50<Ni(i)&Ni(i)<=100plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif100<Ni(i)&Ni(i)<=150plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Ni(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendtext(100,17900,'▲-12.3<Ni(μg/g)<19.9','color','w','backgroundcolor','k');text(100,16900,'★-19.9<=Ni(μg/g)<=50','color','w','backgroundcolor','k');text(100,15900,'●-50<Ni(μg/g)<=100','color','w','backgroundcolor','k');text(100,14900,'◆-100<Ni(μg/g)<150','color','w','backgroundcolor','k');title('Ni的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Pbfig8=figure(8);set(fig8,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高線圖holdonfori=1:X_hangif31<=Pb(i)&Pb(i)<=43plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif43<Pb(i)&Pb(i)<=80plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif80<Pb(i)&Pb(i)<=150plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif150<Pb(i)&Pb(i)<=250plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif250<Pb(i)&Pb(i)<=480plot(x(i),y(i),'s','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Pb(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendendtext(100,17900,'▲-31<Pb(μg/g)<43','color','w','backgroundcolor','k');text(100,16900,'★-43<=Pb(μg/g)<=80','color','w','backgroundcolor','k');text(100,15900,'●-80<Pb(μg/g)<=150','color','w','backgroundcolor','k');text(100,14900,'◆-150<Pb(μg/g)<250','color','w','backgroundcolor','k');text(100,13900,'■-250<Pb(μg/g)<=480','color','w','backgroundcolor','k');title('Pb的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Znfig9=figure(9);set(fig9,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高線圖holdonfori=1:X_hangif69<=Zn(i)&Zn(i)<=97plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif97<Zn(i)&Zn(i)<=160plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif160<Zn(i)&Zn(i)<=500plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif500<Zn(i)&Zn(i)<=1500plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif1