數(shù)學(xué)建模A題 城市表層土壤重金屬污染分析

1、 城市表層土壤重金屬污染分析摘要 對(duì)于問題一，首先采用克里格插值法根據(jù)附件中給的采樣點(diǎn)的濃度數(shù)據(jù)對(duì)城區(qū)內(nèi)每種重金屬濃度值進(jìn)行插值，進(jìn)而繪制每種重金屬的濃度的空間分布圖，對(duì)此為了方便計(jì)算我們借助于Surfer軟件進(jìn)行繪圖，然后我們采用單因子指數(shù)法和內(nèi)梅羅綜合指數(shù)法對(duì)各功能區(qū)的污染程度進(jìn)行評(píng)價(jià)，但是這兩種方法無法從自然異常中分離出人為異常，為了彌補(bǔ)其不足，采用地累積指數(shù)法1對(duì)土壤重金屬污染做進(jìn)一步評(píng)價(jià)，繼而得到不同區(qū)域重金屬的污染程度。評(píng)價(jià)結(jié)果如表9所示。 對(duì)于問題二采用主成分分析法，旨在利用降維的思想，把多指標(biāo)轉(zhuǎn)化為少數(shù)幾個(gè)綜合指標(biāo)， 根據(jù)主成分分析法的一般步驟，首先對(duì)附件中給出的重金屬濃度進(jìn)行

2、標(biāo)準(zhǔn)化處理，然后得到各金屬之間的相關(guān)系數(shù)矩陣，求出相關(guān)系數(shù)矩陣的特征值和向量值，再得到因子的成分矩陣，確定出主成分的個(gè)數(shù)，計(jì)算出各因子的成分得分矩陣，最后通過算綜合主成分中各個(gè)因子的權(quán)重系數(shù)得到污染性較大的因子，最后分析該污染主要原因。結(jié)果如表14和評(píng)價(jià)結(jié)果。 問題三我們用對(duì)流-擴(kuò)散偏微分方程來進(jìn)行描述，對(duì)流擴(kuò)散方程是描述粘性流體運(yùn)動(dòng)的非線性偏微分方程模型，我們將對(duì)流擴(kuò)散方程進(jìn)行簡化，即變?yōu)槎A橢圓形偏微分方程，利用有限插值數(shù)值法估計(jì)出污染物的濃度分布，并與實(shí)際檢測(cè)值相比較，偏差較小的即為污染源的位置及源強(qiáng)，即將問題轉(zhuǎn)化為非線性最優(yōu)化問題，結(jié)果如表所示。 對(duì)于問題四在分析問題三模型的優(yōu)缺點(diǎn)基礎(chǔ)

3、上，為更好的研究城市地質(zhì)環(huán)境的演變模式，將問題三的模型中二維對(duì)流擴(kuò)散偏微分方程擴(kuò)展到三維 ，這樣就可以全面考慮影響模型準(zhǔn)確性的各相關(guān)參數(shù)，通過收集不同的地理，天氣條件下地質(zhì)元素的空間分布信息，進(jìn)而根據(jù)這些信息，建立優(yōu)化模型，即三維方程的參數(shù)，在此基礎(chǔ)上通過模擬仿真進(jìn)行分析。關(guān)鍵字：單因子指數(shù)法 內(nèi)梅羅綜合指數(shù)法 地累積指數(shù)法 主成分分析法 偏微分模型 一 問題重述隨著城市經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷增加，人類活動(dòng)對(duì)城市環(huán)境質(zhì)量的影響日顯突出。對(duì)城市土壤地質(zhì)環(huán)境異常的查證，以及如何應(yīng)用查證獲得的海量數(shù)據(jù)資料開展城市環(huán)境質(zhì)量評(píng)價(jià)，研究人類活動(dòng)影響下城市地質(zhì)環(huán)境的演變模式，日益成為人們關(guān)注的焦點(diǎn)。

4、按照功能劃分，城區(qū)一般可分為生活區(qū)、工業(yè)區(qū)、山區(qū)、主干道路區(qū)及公園綠地區(qū)等，分別記為1類區(qū)、2類區(qū)、5類區(qū)，不同的區(qū)域環(huán)境受人類活動(dòng)影響的程度不同?，F(xiàn)對(duì)某城市城區(qū)土壤地質(zhì)環(huán)境進(jìn)行調(diào)查。為此，將所考察的城區(qū)劃分為間距1公里左右的網(wǎng)格子區(qū)域，按照每平方公里1個(gè)采樣點(diǎn)對(duì)表層土（0-10厘米深度）進(jìn)行取樣、編號(hào)，并用GPS記錄采樣點(diǎn)的位置。應(yīng)用專門儀器測(cè)試分析，獲得了每個(gè)樣本所含的多種化學(xué)元素的濃度數(shù)據(jù)。另一方面，按照2公里的間距在那些遠(yuǎn)離人群及工業(yè)活動(dòng)的自然區(qū)取樣，將其作為該城區(qū)表層土壤中元素的背景值。附件1列出了采樣點(diǎn)的位置、海拔高度及其所屬功能區(qū)等信息，附件2列出了8種主要重金屬元素在采樣點(diǎn)處的

5、濃度，附件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)境的演變模式，還應(yīng)收集什么信息？有了這些信息，如何建立模型解決問題？ 二 問題分析 對(duì)于問題一，要根據(jù)問題附件中給出的采樣點(diǎn)以及在采樣點(diǎn)處8種主要重金屬元素的信息，研究各金屬元素的空間分布，并分析不同功能區(qū)的污染程度，附件中只給出了金屬元素在采樣點(diǎn)處的濃

6、度，要給出在整個(gè)城區(qū)內(nèi)的空間分布，需要通過插值得到更加密集的濃度分布值，由于城區(qū)內(nèi)樣本點(diǎn)間的重金屬濃度存在空間相關(guān)性，因此可以通過克立格法進(jìn)行插值。在這里我們借助Surfer8.0軟件，也是通過克里格插值的方法生成重金屬含量數(shù)據(jù)，并繪制等值線圖。由于還需要分析不同區(qū)域重金屬的污染程度，污染程度需要抽象為污染指數(shù)，通過比較污染指數(shù)與國家背景值，確定污染程度較大的金屬元素，針對(duì)這個(gè)問題我們采用單因子指數(shù)法和內(nèi)梅羅綜合指數(shù)法1，但是這兩種方法無法從自然異常中分離出人為異常，為了彌補(bǔ)其不足，采用地累積指數(shù)法3對(duì)土壤重金屬污染做進(jìn)一步評(píng)價(jià)，繼而得到不同區(qū)域重金屬的污染程度。 對(duì)于問題二要求通過數(shù)據(jù)說明重

7、金屬污染的主要原因，現(xiàn)有資料表明，某些重金屬空間分布具有相關(guān)性，相關(guān)性較大的重金屬可能在成因和來源上有一定的聯(lián)系。在此我們選用主成分分析法來說明。主成分分析也稱主分量分析，旨在利用降維的思想，把多指標(biāo)轉(zhuǎn)化為少數(shù)幾個(gè)綜合指標(biāo)。在實(shí)證問題研究中，為了全面、系統(tǒng)地分析問題，我們必須考慮眾多影響因素。這些涉及的因素一般稱為指標(biāo)，在多元統(tǒng)計(jì)分析中也稱為變量。因?yàn)槊總€(gè)變量都在不同程度上反映了所研究問題的某些信息，并且指標(biāo)之間彼此有一定的相關(guān)性，因而所得的統(tǒng)計(jì)數(shù)據(jù)反映的信息在一定程度上有重疊。在用統(tǒng)計(jì)方法研究多變量問題時(shí)，變量太多會(huì)增加計(jì)算量和增加分析問題的復(fù)雜性，人們希望在進(jìn)行定量分析的過程中，涉及的變量

8、較少，得到的信息量較多。 根據(jù)主成分分析法的一般步驟，首先對(duì)附件中給出的重金屬濃度進(jìn)行標(biāo)準(zhǔn)化處理，然后得到各金屬之間的相關(guān)系數(shù)矩陣，求出相關(guān)系數(shù)矩陣的特征值和向量值，再得到因子的成分矩陣，確定出主成分的個(gè)數(shù)，計(jì)算出各因子的成分得分矩陣，最后通過算綜合主成分中各個(gè)因子的權(quán)重系數(shù)得到污染性較大的因子，最后分析該污染主要原因。 對(duì)于問題三要求分析重金屬污染物的傳播特征，建立模型并分析污染源的位置，對(duì)此我們用對(duì)流-擴(kuò)散偏微分方程來進(jìn)行描述，對(duì)流擴(kuò)散方程是描述粘性流體運(yùn)動(dòng)的非線性偏微分方程模型，可以刻畫很多自然現(xiàn)象，如污染物的擴(kuò)散、降解、流體流動(dòng)與傳熱等，我們將對(duì)流擴(kuò)散方程進(jìn)行簡化，即變?yōu)槎A橢圓形偏微

9、分方程，利用有限插值數(shù)值法估計(jì)出污染物的濃度分布，并與實(shí)際檢測(cè)值相比較，偏差較小的即為污染源的位置及源強(qiáng)，即將問題轉(zhuǎn)化為非線性最優(yōu)化問題。 對(duì)于問題四在分析問題三模型的優(yōu)缺點(diǎn)基礎(chǔ)上，為更好的研究城市地質(zhì)環(huán)境的演變模式，應(yīng)該確定收集什么信息，有了這些信息如何建立模型來解決問題。首先將問題三的模型中二維對(duì)流擴(kuò)散偏微分方程擴(kuò)展到三維 ，這樣就可以全面考慮影響模型準(zhǔn)確性的各相關(guān)參數(shù)，為了得到參數(shù)的準(zhǔn)確表達(dá)式，通過收集不同的地理，天氣條件下地質(zhì)元素的空間分布信息，進(jìn)而根據(jù)這些信息，建立優(yōu)化模型，即三維方程的參數(shù)，在此基礎(chǔ)上通過模擬仿真進(jìn)行分析。 三 問題假設(shè)1.問題一中國家評(píng)價(jià)標(biāo)準(zhǔn)采用二級(jí)標(biāo)準(zhǔn)。2.假設(shè)

10、全部采樣點(diǎn)在自然情況下產(chǎn)生，排除其他因素的干擾。3.假設(shè)該城區(qū)沒有受到災(zāi)難性的污染。4.觀測(cè)值能夠反映一段時(shí)間內(nèi)保持穩(wěn)定的濃度分布。5.假設(shè)污染源是均勻向外擴(kuò)散的。 四 符號(hào)說明 污染物類別 土壤中污染物實(shí)測(cè)濃度 土壤中污染物的背景值 土壤單因子污染指數(shù) 為綜合污染指數(shù) 元素污染指數(shù)平均值 元素污染指數(shù)中的最大值 地累積指數(shù) 各個(gè)區(qū)域元素的平均值 元素地球化學(xué)背景值 各地巖石差異導(dǎo)致的背景值變動(dòng)系數(shù) 第個(gè)采樣點(diǎn)的第種重金屬的濃度值 第種重金屬的均值 第種重金屬的標(biāo)準(zhǔn)差 二元素相關(guān)系數(shù) 五 建模與求解5.1問題一建模與求解 為了得到更加形象具體的重金屬空間分布圖，在這里我們借助Surfer8.0

11、軟件，也是通過克里格插值的方法生成重金屬含量數(shù)據(jù)，并繪制等值線圖。 圖1 城區(qū)地形三維圖 圖2各功能區(qū)地形圖 圖3 城區(qū)濃度分布圖圖4 城區(qū)濃度分布圖 圖5 城區(qū)濃度分布圖 圖6 城區(qū)濃度分布圖 圖7 城區(qū)濃度分布圖 圖8 城區(qū)濃度分布圖 圖9 城區(qū)濃度分布圖 圖10 城區(qū)濃度分布圖5.1.1單因子指數(shù)法 我們對(duì)該城區(qū)不同區(qū)域的污染程度進(jìn)行評(píng)價(jià)，由于每個(gè)區(qū)域都存在八種重金屬，所以我們對(duì)每個(gè)區(qū)域的八種重金屬含量分別進(jìn)行評(píng)價(jià)，評(píng)價(jià)模式分為單因子指數(shù)法和內(nèi)梅羅綜合指數(shù)法。首先建立單因子污染指數(shù)式為： （1）式中為土壤中污染物i的污染分指數(shù)，為土壤中污染物實(shí)測(cè)濃度，為土壤中污染物的背景值。然后通過求解

12、上述單項(xiàng)污染指數(shù)式得到五個(gè)區(qū)域種污染元素的分指數(shù)，如表1（表2、表3、表4、表5）所示為生活區(qū)污染分指數(shù)。 表1 各區(qū)域重金屬實(shí)測(cè)平均濃度編號(hào)As(g/g)Cd(ng/g)Cr(g/g)Cu(g/g)Hg(ng/g)Ni(g/g)Pb(g/g)Zn(g/g)16.27289.9669.0249.4093.0418.3469.11237.0127.25393.1153.41127.54642.3619.8193.04277.9334.04152.3238.9617.3240.9615.4536.5673.2945.71360.0158.0562.21446.8217.6263.53242.855

13、6.26280.5443.6430.19114.9915.2960.71154.24 表2 各區(qū)域單因子污染分指數(shù)編號(hào)As(g/g)Cd(ng/g)Cr(g/g)Cu(g/g)Hg(ng/g)Ni(g/g)Pb(g/g)Zn(g/g)12.012.232.013.742.661.492.233.4322.013.021.729.6618.351.613.004.0331.121.171.261.311.171.261.181.0641.592.771.874.7112.771.432.053.5251.742.161.412.293.291.241.962.24 表3 土壤單因子污染指數(shù)評(píng)價(jià)標(biāo)

14、準(zhǔn)分級(jí)指數(shù)質(zhì)量指數(shù)分級(jí)指數(shù)質(zhì)量指數(shù) 清潔 輕污染潛在污染 重污染5.1.2內(nèi)梅羅綜合指數(shù)法 單項(xiàng)污染指數(shù)得到每個(gè)區(qū)域各種重金屬污染程度，為了得到每個(gè)區(qū)域的綜合污染程度，我們利用內(nèi)梅羅綜合指數(shù)法求解。內(nèi)梅羅綜合指數(shù)式為： （2） 式中：為綜合污染指數(shù), 為各元素污染指數(shù)平均值, 為各元素污染指數(shù)中的最大值。 通過求解內(nèi)梅羅綜合指數(shù)式得到每個(gè)區(qū)域的綜合污染指數(shù)（見表7，表8） 表4 各區(qū)域綜合污染指數(shù)I最大I平均綜合污染指數(shù)P生活區(qū)1.660550.596230.882173工業(yè)區(qū)5.2449541.4343082,718,767山區(qū)0.1889910.0709550.100936主干道區(qū)2.24

15、81650.8782921.206819公園綠地0.8855290.3774380.481306 表5 土壤內(nèi)梅羅綜合污染指數(shù)評(píng)價(jià)標(biāo)準(zhǔn)等級(jí)內(nèi)梅羅綜合污染指數(shù)污染等級(jí)0.7清潔0.71警戒線12輕度污染23重度污染3重污染5.1.3地累積指數(shù)法雖然單因子指數(shù)法和內(nèi)梅羅指數(shù)法均能對(duì)土壤重金屬污染程度進(jìn)行較為全面的評(píng)價(jià)，但無法從自然異常中分離人為異常，然而地累積指數(shù)法彌補(bǔ)了這項(xiàng)不足。地累積指數(shù)法表達(dá)式為： （3）式中：為地累積指數(shù)，是各個(gè)區(qū)域元素的平均值，是元素地球化學(xué)背景值，為各地巖石差異導(dǎo)致的背景值變動(dòng)系數(shù)（一般取值為1.5）。 通過對(duì)地累積指數(shù)式的求解得到每個(gè)區(qū)域各種重金屬的污染程度（見表9）

16、 表6 每個(gè)區(qū)域各種重金屬的污染程度評(píng)價(jià)結(jié)果 經(jīng)過內(nèi)梅羅綜合指數(shù)法和地累積指數(shù)法對(duì)該城不同區(qū)域重金屬的污染程度分析，得到生活區(qū)中、元素屬于輕度污染，元素成為無污染元素，、污染為中度污染。綜合污染等級(jí)為警戒級(jí)。工業(yè)區(qū)中、元素屬于輕度污染，、元素成為中度污染，則屬于強(qiáng)度污染。綜合污染等級(jí)為嚴(yán)重污染。山區(qū)、八種元素基本都對(duì)當(dāng)?shù)丨h(huán)境無污染，綜合污染等級(jí)為清潔級(jí)。主干道路區(qū)、 等元素基本為輕度污染，只有為強(qiáng)度污染，綜合污染等級(jí)為輕度污染。公園綠地區(qū)大部分元素為輕度污染，有很少地方存在中度污染，綜合污染等級(jí)為清潔級(jí)。5.2問題二的建模與求解 根據(jù)問題的分析，某些重金屬空間分布的含量也具有一定的相關(guān)性，相關(guān)

17、性 較大的金屬可能在成因和來源上有一定的聯(lián)系。因此，需要對(duì)問題中所給出的8種重金屬元素之間的關(guān)系進(jìn)行分析。在這里我們選用主成分分析法來求解這個(gè)問題。 主成分分析也稱主分量分析，旨在利用降維的思想，把多指標(biāo)轉(zhuǎn)化為少數(shù)幾個(gè)綜合指標(biāo)。在實(shí)證問題研究中，為了全面、系統(tǒng)地分析問題，我們必須考慮眾多影響因素。這些涉及的因素一般稱為指標(biāo)，在多元統(tǒng)計(jì)分析中也稱為變量。因?yàn)槊總€(gè)變量都在不同程度上反映了所研究問題的某些信息，并且指標(biāo)之間彼此有一定的相關(guān)性，因而所得的統(tǒng)計(jì)數(shù)據(jù)反映的信息在一定程度上有重疊。在用統(tǒng)計(jì)方法研究多變量問題時(shí)，變量太多會(huì)增加計(jì)算量和增加分析問題的復(fù)雜性，人們希望在進(jìn)行定量分析的過程中，涉及的

18、變量較少，得到的信息量較多。5.2.1對(duì)原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化 對(duì)原始數(shù)據(jù)標(biāo)準(zhǔn)化處理的一般公式為： （4）其中為第個(gè)采樣點(diǎn)的第種重金屬的濃度值，分別表示第種重金屬的均值和標(biāo)準(zhǔn)差。5.2.2求重金屬含量相關(guān)系數(shù)矩陣 對(duì)于任意的兩種重金屬元素與元素在第個(gè)樣本中的含量值分別為和（），則二元素相關(guān)系數(shù)的一般計(jì)算公式為 （5） 其中 式中分別為元素與元素樣本點(diǎn)觀測(cè)值的算術(shù)平均值。 首先我們通過將生活區(qū)的重金屬元素原始數(shù)據(jù)導(dǎo)入軟件進(jìn)行主成分分析處理可以得到各重金屬元素之間的相關(guān)系數(shù)如下表所示。（以生活區(qū)為例） 表17 各個(gè)因子間的相關(guān)系數(shù)矩陣5.2.3求相關(guān)系數(shù)矩陣的特征值和特征向量 利用相關(guān)系數(shù)矩陣可求出相

19、應(yīng)的因子特征值和累計(jì)貢獻(xiàn)率，用SPSS統(tǒng)計(jì)軟件計(jì)算可得到生活區(qū)的各因子特征值和成分矩陣結(jié)果如下表所示。 表8 解釋的總方差 表9 成分矩陣（載荷矩陣） 主成分個(gè)數(shù)提取原則為主成分對(duì)應(yīng)的特征值大于1的前m個(gè)主成分，通過表21，可知提取3個(gè)主成分。用表16中數(shù)據(jù)除以主成分對(duì)應(yīng)特征值便得到三個(gè)成分中每個(gè)指標(biāo)對(duì)應(yīng)的系數(shù)，即特征向量A1，A2，A3，將特征向量與標(biāo)準(zhǔn)化后數(shù)據(jù)相乘得到主成分表達(dá)式：F1=A1*ZX1; F2=A2*ZX2; F3=A3*ZX3。為了方便計(jì)算我們?cè)谶@里用SPSS解決這個(gè)問題得到成分得分系數(shù)矩陣。 表10 各因子得分系數(shù)矩陣以每個(gè)主成分所對(duì)應(yīng)的特征值占所提取主成分總的特征值之和

20、的比例作為權(quán)重計(jì)算主成分綜合模型，即用第一主成分F1中每個(gè)指標(biāo)所對(duì)應(yīng)的系數(shù)乘上第一主成分F1所對(duì)應(yīng)的貢獻(xiàn)率再除以所提取的三個(gè)主成分貢獻(xiàn)率之和，然后加上第二主成分F2對(duì)應(yīng)系數(shù)乘以其貢獻(xiàn)率百分比，再加上第三主成分F3對(duì)應(yīng)系數(shù)乘以其貢獻(xiàn)率百分比，就得到了綜合主成分F中每個(gè)指標(biāo)所對(duì)應(yīng)的權(quán)重系數(shù)A。根據(jù)權(quán)重系數(shù)A值大小降序排列，系數(shù)大的對(duì)應(yīng)的重金屬元素就是主要污染的主要原因。這里只將生活區(qū)列舉出來，其它區(qū)域分析方法與上述方法相同，就不做過多敘述，其他區(qū)域表格見附錄。 表11 各區(qū)域因子的權(quán)重系數(shù)生活區(qū)權(quán)重系數(shù)工業(yè)區(qū)權(quán)重系數(shù)山區(qū)權(quán)重系數(shù)Cr0.34Ni0.34Cu0.35Zn0.34Zn0.33As0.2

21、3Ni0.27As0.31Cr0.23Pb0.22Pb0.31Ni0.22Cd0.21Cd0.29Hg0.21Cu0.2Cr0.29Zn0.17Hg0.11Cu0.23Pb0.06As0.1Hg0.22Cd0.04主干道區(qū)權(quán)重系數(shù)公園綠地權(quán)重系數(shù)Pb0.38Pb0.36Cd0.34Cu0.27Cu0.29Hg0.27Ni0.27Zn0.26Zn0.27Cd0.23Cr0.26Cr0.16Hg0.23As0.14As0.06Ni0.08結(jié)果評(píng)價(jià) 對(duì)該城區(qū)不同區(qū)域采集的土壤、元素進(jìn)行主成分分析，通過主成分分析法得到: 1> 生活區(qū)元素、在該城市土壤中的含量明顯高于國家背景值。這是由于生活區(qū)存

22、在大量的公路，公路兩側(cè)含鉛汽油的燃燒和汽車輪胎的磨損的粉塵會(huì)增加土壤中、等元素的含量，而且由于生活區(qū)人類活動(dòng)比較頻繁，日常生活用品丟棄后成為垃圾也會(huì)導(dǎo)致這三種元素的增加。2> 工業(yè)區(qū)主要污染元素為、。因?yàn)楣I(yè)區(qū)“三廢”排放，采礦和冶煉會(huì)增加這些元素在土壤中的含量。 3> 山區(qū)全部元素幾乎都在國家背景值以下，只有元素為主要污染元素。這是由于山區(qū)遠(yuǎn)離城市和工業(yè)區(qū)的污染源，只有一部分金屬礦山的開采會(huì)導(dǎo)致一些金屬元素的外露。 4> 主干道路區(qū)、元素都高于背景值，其余元素都低于背景值。機(jī)動(dòng)車尾氣排放既是城市大氣的主要污染源，也顯著引起公路兩側(cè)土壤的重金屬污染，汽車汽油、發(fā)動(dòng)機(jī)、輪胎、潤

23、滑油和鍍金部分都能燃燒或磨損而釋放出、元素。5> 公園綠地區(qū)中為主要重金屬元素，、為次要重金屬元素。公園綠地區(qū)農(nóng)藥化肥的使用會(huì)不同程度的影響重金屬的污染，而且地下水灌溉和塑料薄膜會(huì)增加、元素的土壤含量。不難發(fā)現(xiàn)每個(gè)區(qū)域都受到不同程度的元素的污染，這是因?yàn)殒k元素主要以硫化鎘形式儲(chǔ)存于鋅礦、鉛鋅礦和銅鉛鋅礦中，土壤鎘主要來源于鋅礦、鉛礦的冶煉、合金、電鍍、化工廠等排放的廢水，工業(yè)固廢堆放，含有廢舊電池的生活垃圾滲濾，污泥施肥以及過量或不恰當(dāng)?shù)厥褂没兽r(nóng)藥等。總體來說：工業(yè)化程度越高的地區(qū)污染越嚴(yán)重，市區(qū)高于遠(yuǎn)郊和農(nóng)村，地表高于地下，污染區(qū)污染時(shí)間越長重金屬積累就越多，以大氣傳播媒介土壤重金屬

24、污染土壤具有很強(qiáng)的疊加性,熟化程度越高重金屬含量就越高。5.3 問題三模型建立與求解 通過對(duì)問題二的分析，我們知道了重金屬的傳播途徑，比如重金屬和的污染主要為汽車尾氣排放所致，的污染為燃煤引起，它們都是經(jīng)過空氣來傳播的，為此用對(duì)流-擴(kuò)散偏微分方程來描述。主要是由于工業(yè)廢水排放污染，經(jīng)過地表徑流的傳播，則可以應(yīng)用對(duì)流-擴(kuò)散偏微分方程來進(jìn)行描述。對(duì)流-擴(kuò)散方程是描述粘性流體運(yùn)動(dòng)的非線性方程的線性化偏微分方程模型，能夠刻畫很多自然現(xiàn)象，如污染物的擴(kuò)散，降解，流體流動(dòng)與傳熱和電化學(xué)反應(yīng)等。 首先建立對(duì)流-擴(kuò)散偏微分方程傳播模型。然后，依據(jù)問題所給的條件，可以知道這是對(duì)流-擴(kuò)散偏微分方程的反問題。在求解

25、時(shí)，將對(duì)流-擴(kuò)散偏微分方程進(jìn)行簡化變?yōu)槎A橢圓型偏微分方程。利用橢圓型偏微分方程的有限差分?jǐn)?shù)值解法正向估計(jì)出污染物的濃度分布，再代入不同的位置和強(qiáng)度參數(shù)，與實(shí)際檢測(cè)值相比較，偏差最小的即為污染源的位置及源強(qiáng)，即將問題轉(zhuǎn)化為非線性最優(yōu)化問題，利用中工具箱求解橢圓型微分方程，則可得到各重金屬元素的污染源位置及源強(qiáng)。 對(duì)流-擴(kuò)散偏微分方程模型的建立一般的對(duì)流-擴(kuò)散偏微分方程模型為 （6）其中為有界區(qū)域，為邊界上的外法向量，為軸方向的流速，為軸方向的流速，為擴(kuò)散系數(shù)，為污染物的降解率，為污染物的濃度分布，和分別表示多個(gè)污染點(diǎn)和污染強(qiáng)度，為狄拉克函數(shù)，即單位脈沖函數(shù)。 5.3.2對(duì)流-擴(kuò)散偏微分方程反問

26、題的求解若已知的分布，并且已知和那么源項(xiàng)識(shí)別的反問題就是根據(jù)這些已知的分布確定源項(xiàng)，即確定污染源的位置和污染強(qiáng)度。根據(jù)模型的假設(shè)，觀測(cè)值為一段時(shí)間內(nèi)穩(wěn)定的濃度分布，即二階穩(wěn)態(tài)對(duì)流擴(kuò)散方程。根據(jù)相關(guān)文獻(xiàn)可知擴(kuò)散系數(shù)=，污染物降解率=0.001，暫時(shí)不考慮橫縱軸的流速，使軸方向的流速=0，軸方向的流速=0.于是則可將上述的偏微分方程簡化為 （7）該方程為二階橢圓型偏微分方程。要在污染區(qū)域中大概找出污染范圍，首先選擇污染源初始位置，以及污染強(qiáng)度。利用橢圓型偏微分方程的有限差分?jǐn)?shù)值解正向估計(jì)出污染物的濃度分布，。 的分布與初始位置，以及污染強(qiáng)度的取值有關(guān)。用函數(shù)表示在初始位置，以及污染強(qiáng)度的條件下測(cè)得

27、的污染物濃度分布，為第次估計(jì)值。將所得到的污染物濃度與實(shí)際檢測(cè)值相比較，差距最小的即為最接近污染源的位置。即將問題轉(zhuǎn)換為非線性最優(yōu)化問題： （8）這里為選取的樣點(diǎn)個(gè)數(shù)。 適時(shí)調(diào)整參數(shù)，使上面的問題達(dá)到最小值。使其達(dá)到最小值的，即為污染源的位置坐標(biāo)和污染強(qiáng)度。 5.3.3對(duì)流-擴(kuò)散方程解的仿真模擬對(duì)于簡化的對(duì)流-擴(kuò)散方程化為標(biāo)準(zhǔn)的二階橢圓型偏微分方程： （9）用函數(shù)來近似表示單位脈沖函數(shù)，從而可得到偏微分方程為： 通過以上的分析，結(jié)合問題（1）的模型所得到的重金屬污染分布結(jié)果可知，污染濃度最高的區(qū)域最有可能是污染源。因此搜索污染源時(shí)以差值所得到的濃度最高點(diǎn)作為搜索的初始值，并在該值附近進(jìn)行搜索。

28、對(duì)每個(gè)可能的污染點(diǎn)使用中的工具來求解橢圓型偏微分方程的解。比較不同污染點(diǎn)和源強(qiáng)下的污染分布和實(shí)際檢測(cè)的污染濃度，選擇兩者誤差最小的污染源位置和源強(qiáng)。用此方法可求得八種重金屬的污染源位置和污染濃度如下： 表12 各種金屬污染源坐標(biāo)和濃度值重金屬橫坐標(biāo)/縱坐標(biāo)/濃度值(g/g)As181341004624.544913289251519.27Cd 236649790725.993921242486542.6413Cr 32996018320.3774Cu 23833692447.5421Hg 270822951672.21410023671508.11518591211326.8Ni 329960

29、1832.9542Pb 47774897165.78Zn 1324470561289.56913245821156.75.4問題四模型建立與求解 在解決問題三的模型中，將復(fù)雜的對(duì)流-擴(kuò)散方程簡化得到污染物傳播的基本模型，便于直觀地從解空間中看到傳播趨勢(shì)，從而便于反問題求解污染源，但是沒有考慮橫縱軸的流速，對(duì)于隨位置變化的擴(kuò)散系數(shù)，以及降解率等以常數(shù)值來描述，具有局限性。無法滿足對(duì)城市地質(zhì)環(huán)境研究的需求。需要進(jìn)一步改進(jìn)。5.4.1擴(kuò)展二維對(duì)流擴(kuò)散方程到三維 首先要研究哪些因素對(duì)地質(zhì)環(huán)境的演變影響比較大，比如不僅受地域地形的影響，還受到天氣(如風(fēng))、雨水沖刷和溫度的影響，為了更準(zhǔn)確的研究城市地質(zhì)環(huán)

30、境的演變，還需要收集發(fā)地區(qū)在不同天氣條件下各種重金屬的濃度分布。如在同一區(qū)域內(nèi)雨天前的濃度分布與雨天后的濃度分布情況，在不同季節(jié)、不同溫度條件下的濃度分布，以及風(fēng)的因素隨污染物濃度的影響。 因此，要實(shí)現(xiàn)這個(gè)目的，就需要將問題三中的二維對(duì)流-擴(kuò)散模型擴(kuò)展到三維的情形，即有三維污染物的擴(kuò)散模型為 （10）其中為污染物濃度；為軸方向的流速；為軸方向的流速；為軸方向的流速，表示擴(kuò)散系數(shù)；為自凈化系數(shù)，表示沉降、稀釋、化學(xué)及生物反應(yīng)等，如當(dāng)污染物為重金屬時(shí)表示重金屬形態(tài)轉(zhuǎn)化系數(shù)，當(dāng)污染物為有機(jī)物時(shí)表示有機(jī)物降解系數(shù)，亦可是各因素的綜合；表示污染物濃度隨時(shí)間的變化項(xiàng)。 在均勻穩(wěn)定流中，三維對(duì)流-擴(kuò)散模型的

31、解為 （11）其中為源強(qiáng)。5.4.2 三維對(duì)流-擴(kuò)散偏微分方程參數(shù)的確定 在這里主要是估計(jì)出更準(zhǔn)確的分布參數(shù)，從而確定污染物的分布個(gè)數(shù)，并能夠有效地模擬實(shí)際中地質(zhì)模型的演變模式，問題就轉(zhuǎn)換為準(zhǔn)確估計(jì)三維對(duì)流-擴(kuò)散方程的參數(shù)，進(jìn)而用以模擬城市的地質(zhì)演變模式。 在坡度和溫度基本相同的條件下，考慮風(fēng)的影響。某些重金屬，如和主要來自于燃煤或工業(yè)廢氣引起的徑流污染，這就使得風(fēng)向和風(fēng)速對(duì)這些重金屬分布的影響很大，在大氣中重金屬會(huì)隨著風(fēng)的方向而漂移，從而改變了一定區(qū)域內(nèi)的重金屬含量濃度。 在不同的風(fēng)向量條件下，收集區(qū)域土地的分布，這里為風(fēng)向量，即有大小，又有方向。軸方向的流速為 （12）將上式代入到（15.

32、4）式中得到的數(shù)學(xué)模型為 當(dāng)函數(shù)表示不同形式時(shí)，利用解析解得到不同的重金屬濃度分布。將由不同函數(shù)所得到的污染物濃度物空間濃度分布與觀測(cè)值進(jìn)行比較，則可將問題轉(zhuǎn)化為求解優(yōu)化問題，即利用多項(xiàng)式逼近函數(shù)，從而確定風(fēng)對(duì)參數(shù)影響的近似函數(shù)表達(dá)式。 首先考慮風(fēng)速對(duì)軸方向的影響，對(duì)模型進(jìn)行仿真模擬得到的污染物濃度的空間分布趨勢(shì)，說明風(fēng)對(duì)污染物濃度的影響。 同理對(duì)于軸和軸方向，以及地形和溫度等因素對(duì)三維對(duì)流-擴(kuò)散方程的影響，類似的也可以用上面的模型進(jìn)行估計(jì)修正，最后綜合不同因素對(duì)參數(shù)的影響則可以確定參數(shù)形式的三維對(duì)流-擴(kuò)散模型。從而可以從該模型出發(fā)，更好的研究城市地質(zhì)演變的模式，有針對(duì)性的對(duì)城市的地質(zhì)環(huán)境問題

33、進(jìn)行處理，這樣會(huì)具有很高的實(shí)際應(yīng)用價(jià)值。 六 模型評(píng)價(jià) 6.1評(píng)價(jià) 對(duì)于問題一,利用重金屬濃度在空間的相關(guān)性，采樣普通克里格插值法，繪出各種重金屬的濃度空間分布，然后利用地質(zhì)積累法和修正的內(nèi)梅羅綜合指數(shù)法對(duì)各功能區(qū)的污染程度進(jìn)行分析，給出了各功能區(qū)被重金屬污染的情況，同時(shí)給出了每種重金屬的污染指數(shù)分布，以及在各功能區(qū)的污染指數(shù)直方圖，多方面分析了重金屬的污染程度，但克里格插值法有一定的局限性，在很大的程度上受局部采樣點(diǎn)含量的影響較大，故對(duì)采樣點(diǎn)所具有的代表性要求較高。 對(duì)于問題二，采用主成分分析法，用較少的有代表性的因子來說明眾多因素變量所提取的主要信息，提取出了6個(gè)因子后計(jì)算其得分，這個(gè)模型

34、在一定程度上降低了計(jì)算的復(fù)雜度，并用合理的方法對(duì)污染源進(jìn)行了分析，其結(jié)果與實(shí)際情況基本相符，由于問題中所提供的數(shù)據(jù)是以1公里為網(wǎng)格進(jìn)行取樣的，沒有給出具體的地理分布，所以在分析污染源的時(shí)是在插值的基礎(chǔ)上進(jìn)行分析的，有一定的模糊性，一些污染物的主要來源是與實(shí)際情況相結(jié)合確定出來的。 對(duì)于問題三，建立了二維對(duì)流擴(kuò)散偏微分方程來描述重金屬污染的傳播特征，這是比較合理的，在求解模型時(shí)將其簡化為橢圓形偏微分方程模型，并將問題轉(zhuǎn)化為非線性最優(yōu)解，則較容易地得到問題的結(jié)果但在簡化過程中忽略了一些變量因素對(duì)污染傳播的影響，具有一定的局限性。 對(duì)于問題四，將二維對(duì)流擴(kuò)散方程擴(kuò)展到三維，并對(duì)各種變化參數(shù)進(jìn)行了討論

35、，使更充分地描述各參數(shù)對(duì)城市地質(zhì)演變過程的影響，但沒有給出更詳細(xì)的在各種變化參數(shù)影響下的傳播仿真模擬結(jié)果。6.2模型的改進(jìn)方向 在問題三的模型中，在應(yīng)用對(duì)流-擴(kuò)散方程的時(shí)候，相關(guān)參數(shù)的確定是根據(jù)專家咨詢法給出的，有一定的主觀性，在實(shí)際中為了更準(zhǔn)確地模擬重金屬的污染擴(kuò)散，應(yīng)該根據(jù)一定的科學(xué)方法適當(dāng)?shù)匦拚齾?shù)。對(duì)于偏微分方程的求解，可以應(yīng)用更準(zhǔn)確的算法進(jìn)行求解。在實(shí)際問題中，對(duì)于偏微分方程的求解可以應(yīng)用蒙特卡洛，等方法，在大量數(shù)據(jù)的基礎(chǔ)上對(duì)偏微分方程的參數(shù)進(jìn)行反演，這樣可以提高精度，從而更好地應(yīng)用于研究實(shí)際中城市地質(zhì)環(huán)境模式。 七 參考文獻(xiàn)1何瑤，劉穎，陳玲，上海市表層土壤中金屬元素含量及其分布特

36、征 J.中 國可持續(xù)發(fā)展論壇?？?008。2 趙靜、但琦，數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)M,北京：高等教育出版，2008。3 李尚志，數(shù)學(xué)建模競(jìng)賽教程M，南京：江蘇教育出版社，1996。4 王樹禾，數(shù)學(xué)模型基礎(chǔ)M,臺(tái)肥：中國科學(xué)技術(shù)大學(xué)出版社，1996。5 姜啟源，數(shù)學(xué)模型M，北京：高等教育出版社，1987。6 薛定宇,高等應(yīng)用數(shù)學(xué)問題的MATLAB求解M，北京：清華大學(xué)出版社，2004。7 陳杰，MATLAB寶典M，北京：電子工業(yè)出版社，2007年。 八 附錄程序一：load('chenshiturangjinshu.mat');%x,y=meshgrid(x_zuobiao,y_zu

37、obiao);x=x_zuobiao;y=y_zuobiao;z=haiba;X,Y=meshgrid(0:50:30000,0:50:20000);Z=griddata(x,y,z,X,Y);meshz(X,Y,Z);shading interp;程序二：function pdemodelpde_fig,ax=pdeinit;pdetool('appl_cb',1);set(ax,'DataAspectRatio',15000 10000 1);set(ax,'PlotBoxAspectRatio',1 1 1);set(ax,'XLi

38、m',0 30000);set(ax,'YLim',0 20000);set(ax,'XTickMode','auto');set(ax,'YTickMode','auto');pdetool('gridon','on');% Geometry description:pderect(15161.538461538461 19776.923076923071 10207.692307692307 6515.3846153846134,'R1');set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','R1')% Boundary conditions:pdetool('changemode',0)pdesetbd(4,.'dir',.1,.'1',.'0')pdesetbd(3,.'dir',.1,.'1',.