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

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

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

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

4、按照功能劃分，城區(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個采樣點對表層土（0-10厘米深度）進(jìn)行取樣、編號，并用GPS記錄采樣點的位置。應(yīng)用專門儀器測試分析，獲得了每個樣本所含的多種化學(xué)元素的濃度數(shù)據(jù)。另一方面，按照2公里的間距在那些遠(yuǎn)離人群及工業(yè)活動的自然區(qū)取樣，將其作為該城區(qū)表層土壤中元素的背景值。附件1列出了采樣點的位置、海拔高度及其所屬功能區(qū)等信息，附件2列出了8種主要重金屬元素在采樣點處的

5、濃度，附件3列出了8種主要重金屬元素的背景值?，F(xiàn)要求你們通過數(shù)學(xué)建模來完成以下任務(wù)：(1) 給出8種主要重金屬元素在該城區(qū)的空間分布，并分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。(2) 通過數(shù)據(jù)分析，說明重金屬污染的主要原因。(3) 分析重金屬污染物的傳播特征，由此建立模型，確定污染源的位置。(4) 分析你所建立模型的優(yōu)缺點，為更好地研究城市地質(zhì)環(huán)境的演變模式，還應(yīng)收集什么信息？有了這些信息，如何建立模型解決問題？ 二 問題分析 對于問題一，要根據(jù)問題附件中給出的采樣點以及在采樣點處8種主要重金屬元素的信息，研究各金屬元素的空間分布，并分析不同功能區(qū)的污染程度，附件中只給出了金屬元素在采樣點處的濃

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

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

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

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

10、全部采樣點在自然情況下產(chǎn)生，排除其他因素的干擾。3.假設(shè)該城區(qū)沒有受到災(zāi)難性的污染。4.觀測值能夠反映一段時間內(nèi)保持穩(wěn)定的濃度分布。5.假設(shè)污染源是均勻向外擴散的。 四 符號說明 污染物類別 土壤中污染物實測濃度 土壤中污染物的背景值 土壤單因子污染指數(shù) 為綜合污染指數(shù) 元素污染指數(shù)平均值 元素污染指數(shù)中的最大值 地累積指數(shù) 各個區(qū)域元素的平均值 元素地球化學(xué)背景值 各地巖石差異導(dǎo)致的背景值變動系數(shù) 第個采樣點的第種重金屬的濃度值 第種重金屬的均值 第種重金屬的標(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ù)法 我們對該城區(qū)不同區(qū)域的污染程度進(jìn)行評價，由于每個區(qū)域都存在八種重金屬，所以我們對每個區(qū)域的八種重金屬含量分別進(jìn)行評價，評價模式分為單因子指數(shù)法和內(nèi)梅羅綜合指數(shù)法。首先建立單因子污染指數(shù)式為： （1）式中為土壤中污染物i的污染分指數(shù)，為土壤中污染物實測濃度，為土壤中污染物的背景值。然后通過求解

12、上述單項污染指數(shù)式得到五個區(qū)域種污染元素的分指數(shù)，如表1（表2、表3、表4、表5）所示為生活區(qū)污染分指數(shù)。 表1 各區(qū)域重金屬實測平均濃度編號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ù)編號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ù)評價標(biāo)

14、準(zhǔn)分級指數(shù)質(zhì)量指數(shù)分級指數(shù)質(zhì)量指數(shù) 清潔 輕污染潛在污染 重污染5.1.2內(nèi)梅羅綜合指數(shù)法 單項污染指數(shù)得到每個區(qū)域各種重金屬污染程度，為了得到每個區(qū)域的綜合污染程度，我們利用內(nèi)梅羅綜合指數(shù)法求解。內(nèi)梅羅綜合指數(shù)式為： （2） 式中：為綜合污染指數(shù), 為各元素污染指數(shù)平均值, 為各元素污染指數(shù)中的最大值。 通過求解內(nèi)梅羅綜合指數(shù)式得到每個區(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ù)評價標(biāo)準(zhǔn)等級內(nèi)梅羅綜合污染指數(shù)污染等級0.7清潔0.71警戒線12輕度污染23重度污染3重污染5.1.3地累積指數(shù)法雖然單因子指數(shù)法和內(nèi)梅羅指數(shù)法均能對土壤重金屬污染程度進(jìn)行較為全面的評價，但無法從自然異常中分離人為異常，然而地累積指數(shù)法彌補了這項不足。地累積指數(shù)法表達(dá)式為： （3）式中：為地累積指數(shù)，是各個區(qū)域元素的平均值，是元素地球化學(xué)背景值，為各地巖石差異導(dǎo)致的背景值變動系數(shù)（一般取值為1.5）。 通過對地累積指數(shù)式的求解得到每個區(qū)域各種重金屬的污染程度（見表9）

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

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

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

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

20、的比例作為權(quán)重計算主成分綜合模型，即用第一主成分F1中每個指標(biāo)所對應(yīng)的系數(shù)乘上第一主成分F1所對應(yīng)的貢獻(xiàn)率再除以所提取的三個主成分貢獻(xiàn)率之和，然后加上第二主成分F2對應(yīng)系數(shù)乘以其貢獻(xiàn)率百分比，再加上第三主成分F3對應(yīng)系數(shù)乘以其貢獻(xiàn)率百分比，就得到了綜合主成分F中每個指標(biāo)所對應(yīng)的權(quán)重系數(shù)A。根據(jù)權(quán)重系數(shù)A值大小降序排列，系數(shù)大的對應(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é)果評價 對該城區(qū)不同區(qū)域采集的土壤、元素進(jìn)行主成分分析，通過主成分分析法得到: 1> 生活區(qū)元素、在該城市土壤中的含量明顯高于國家背景值。這是由于生活區(qū)存

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

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

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

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

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

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

28、對每個可能的污染點使用中的工具來求解橢圓型偏微分方程的解。比較不同污染點和源強下的污染分布和實際檢測的污染濃度，選擇兩者誤差最小的污染源位置和源強。用此方法可求得八種重金屬的污染源位置和污染濃度如下： 表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ù)雜的對流-擴散方程簡化得到污染物傳播的基本模型，便于直觀地從解空間中看到傳播趨勢，從而便于反問題求解污染源，但是沒有考慮橫縱軸的流速，對于隨位置變化的擴散系數(shù)，以及降解率等以常數(shù)值來描述，具有局限性。無法滿足對城市地質(zhì)環(huán)境研究的需求。需要進(jìn)一步改進(jìn)。5.4.1擴展二維對流擴散方程到三維 首先要研究哪些因素對地質(zhì)環(huán)境的演變影響比較大，比如不僅受地域地形的影響，還受到天氣(如風(fēng))、雨水沖刷和溫度的影響，為了更準(zhǔn)確的研究城市地質(zhì)環(huán)

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

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

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

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

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

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

36、征 J.中 國可持續(xù)發(fā)展論壇專刊，2008。2 趙靜、但琦，數(shù)學(xué)建模與數(shù)學(xué)實驗M,北京：高等教育出版，2008。3 李尚志，數(shù)學(xué)建模競賽教程M，南京：江蘇教育出版社，1996。4 王樹禾，數(shù)學(xué)模型基礎(chǔ)M,臺肥：中國科學(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',.