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

2、報名號的話）： 20001032 所屬學(xué)校（請?zhí)顚懲暾娜?中國人民解放軍國防科學(xué)技術(shù)大學(xué) 參賽隊員 (打印并簽名) ：1. 李騰驥 2. 袁雪強(qiáng) 3. 趙泉樸 指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人 (打印并簽名)： 楊文強(qiáng) 日期： 2011 年 9 月 9 日賽區(qū)評閱編號（由賽區(qū)組委會評閱前進(jìn)行編號）：2011高教社杯全國大學(xué)生數(shù)學(xué)建模競賽編 號 專 用 頁賽區(qū)評閱編號（由賽區(qū)組委會評閱前進(jìn)行編號）：賽區(qū)評閱記錄（可供賽區(qū)評閱時使用）：評閱人評分備注全國統(tǒng)一編號（由賽區(qū)組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進(jìn)行編號）：基于MCMC的城市表層土壤重金屬污染分析摘 要：本文針對某城

3、市城區(qū)表層土壤重金屬污染分析的問題，從各分區(qū)重金屬元素污染的空間分布、污染程度、污染原因及傳播途徑等多方面進(jìn)行了全面的統(tǒng)計與研究，同時采用貝葉斯推理和馬爾科夫鏈蒙特卡洛法（MCMC），利用MATLAB進(jìn)行編程，對所設(shè)計的模型進(jìn)行求解，得到了各重金屬污染源的地理位置坐標(biāo)，最后對城市地質(zhì)環(huán)境演變模式的問題進(jìn)行了探究，得出了初步的結(jié)論。 首先，根據(jù)所給數(shù)據(jù)信息，利用MATLAB及EXCEL作圖，得到了該地區(qū)的地勢圖及各重金屬污染程度的分布圖。 對于問題1，通過統(tǒng)計各重金屬的重污染點在各區(qū)域的分布和計算不同區(qū)域各中金屬污染程度之間的相關(guān)系數(shù)，得到了各重金屬的大致空間分布，并采用改進(jìn)的內(nèi)梅羅指數(shù)法，得出

4、不同區(qū)域重金屬的污染程度。 對于問題2，在問題1所得數(shù)據(jù)及結(jié)論的基礎(chǔ)上，結(jié)合實際情況，對各重金屬污染的主要原因進(jìn)行了較為全面的分析。 對于問題3，通過對污染的傳播特征的研究，基于城市的地質(zhì)特性，得到空氣傳播為污染主要傳播方式的結(jié)論，并由此建立關(guān)于氣體傳播的概率模型，采用貝葉斯推理和馬爾科夫鏈蒙特卡洛法，利用MATLAB進(jìn)行50000次抽樣，得到收斂的坐標(biāo)值，從而計算出各重金屬污染源的地理位置坐標(biāo)分別為：As（2742,7294）、Cd（2382,3693）、Cr（3470,2309）、Cu（2707,2295）、Hg(2708,2294)、Ni（1647,2729）、Pb（2882,3617）

5、、Zn（4152,2299）。 對于問題4，以題中城市土壤重金屬污染的演變規(guī)律為例進(jìn)行了一定深度的探究，并對所需信息及模型設(shè)計方案進(jìn)行了初步的討論。關(guān)鍵詞：內(nèi)梅羅指數(shù)法 相關(guān)系數(shù) 貝葉斯推理 MCMC法 MH算法 一、問題的提出與分析（一）背景介紹隨著城市經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷增加，人類活動對城市環(huán)境質(zhì)量的影響日顯突出。對城市土壤地質(zhì)環(huán)境異常的查證，以及如何應(yīng)用查證獲得的海量數(shù)據(jù)資料開展城市環(huán)境質(zhì)量評價，研究人類活動影響下城市地質(zhì)環(huán)境的演變模式，日益成為人們關(guān)注的焦點。按照功能劃分，城區(qū)一般可分為生活區(qū)、工業(yè)區(qū)、山區(qū)、主干道路區(qū)及公園綠地區(qū)等，分別記為1類區(qū)、2類區(qū)、5類區(qū)，不同的區(qū)域

6、環(huán)境受人類活動影響的程度不同?，F(xiàn)對某城市城區(qū)土壤地質(zhì)環(huán)境進(jìn)行調(diào)查。為此，將所考察的城區(qū)劃分為間距1公里左右的網(wǎng)格子區(qū)域，按照每平方公里1個采樣點對表層土（010 厘米深度）進(jìn)行取樣、編號，并用GPS記錄采樣點的位置。應(yīng)用專門儀器測試分析，獲得了每個樣本所含的多種化學(xué)元素的濃度數(shù)據(jù)。另一方面，按照2公里的間距在那些遠(yuǎn)離人群及工業(yè)活動的自然區(qū)取樣，將其作為該城區(qū)表層土壤中元素的背景值。（二）要解決的問題現(xiàn)給出采樣點的位置、海拔高度及其所屬功能區(qū)等信息，8種主要重金屬元素在采樣點處的濃度及8種主要重金屬元素的背景值，解答以下問題：1. 給出8種主要重金屬元素在該城區(qū)的空間分布，并分析該城區(qū)內(nèi)不同區(qū)域

7、重金屬的污染程度。2. 通過數(shù)據(jù)分析，說明重金屬污染的主要原因。3. 分析重金屬污染物的傳播特征，由此建立模型，確定污染源的位置。4. 分析所建立模型的優(yōu)缺點，并說明為更好地研究城市地質(zhì)環(huán)境的演變模式，還應(yīng)收集什么信息，有了這些信息，如何建立模型解決問題。（三）問題分析整個問題可以分為三類分別進(jìn)行求解：問題1和問題2可將其歸結(jié)為對所給數(shù)據(jù)的統(tǒng)計與分析，可以通過作圖或列表等直觀的形式針對不同金屬和區(qū)域進(jìn)行相關(guān)統(tǒng)計，再對統(tǒng)計的結(jié)果加以比對與分析，即可得出結(jié)論。問題3確定污染源的位置實質(zhì)就是要找到使周圍污染達(dá)到最大的點，這樣可以轉(zhuǎn)化為概率模型，在統(tǒng)計數(shù)據(jù)的基礎(chǔ)上，運(yùn)用貝葉斯推理馬爾科夫鏈蒙特卡洛法可

8、加以計算和確定。問題4開放類問題，可根據(jù)現(xiàn)有資料進(jìn)行回答。二、模型假設(shè)1. 題目中所給出的樣本均具有普遍性及代表性，而非所測地區(qū)中的特殊情況。2. 該地區(qū)的空氣流動沒有顯著的方向性特征。3. 污染物濃度在取樣前后較長一段時間內(nèi)不會有顯著的變化。三、符號說明: 土壤污染元素綜合污染指數(shù): 土壤中污染元素i的實測值: 土壤中污染元素i的評價標(biāo)準(zhǔn)值: 任一閉合曲面: 內(nèi)指定所圍的區(qū)域: 通過S流入的質(zhì)量: 的質(zhì)量: 內(nèi)由于濃度的變化引起的質(zhì)量增加: 衰減系數(shù): 貝葉斯公式模型參數(shù)m: 模型參數(shù)的個數(shù)y: 貝葉斯公式觀測數(shù)據(jù)：馬爾可夫蒙特鏈卡洛法中的分布：接受概率: 目標(biāo)概率密度函數(shù)s : 距離污染源

9、的距離p： 沉降系數(shù)h: 海拔高度：臨界海拔： 重金屬元素吸附系數(shù)四、模型建立與求解（一）空間分布的特征根據(jù)所給數(shù)據(jù)，對采樣點按區(qū)域進(jìn)行重新的排序，由此作出各采樣點的地勢分布圖（圖一），又由各種金屬在各區(qū)域的濃度值，作出金屬污染程度在各區(qū)域的分布圖（以Pb（鉛）為例，見圖二，其它金屬元素分布圖見附件），由此分析可得到各金屬元素在該城區(qū)的空間分布。圖一 采樣點地勢分布圖海拔（m）圖二 Pb污染程度分布圖以Pb為例進(jìn)行分析，經(jīng)統(tǒng)計，Pb的重度污染區(qū)域主要集中于圖的左下方，且Pb的重度污染點在工業(yè)區(qū)占62.8%，另有20.9%重度污染點則分布在左下方的交通區(qū)中，由此可得Pb的污染區(qū)域主要分布于工業(yè)區(qū)

10、和左下方的部分交通區(qū)，對于其它金屬元素，可以用相同的方法，還可以通過計算不同區(qū)域各金屬元素污染程度的相關(guān)系數(shù)，找出各金屬元素污染之間的聯(lián)系，這樣可以有效提高分析效率和準(zhǔn)確性。例如對表一中工業(yè)區(qū)的各金屬元素污染程度相關(guān)系數(shù)（其它區(qū)域見附件）的比較可以得出，Cr與Cu、Cr與Hg、Cu與Hg、Cd與Pb、Cd與Zn、Pb與Zn相關(guān)系數(shù)較大，可以認(rèn)為它們具有相同或相近的空間分布。表一 工業(yè)區(qū)各金屬元素污染程度相關(guān)系數(shù)CdCrCuHgNiPbZnAs0.3286010.3796280.1529320.1812950.6897250.3946630.517731Cd-0.5409850.5664580.

11、5331710.4886660.8292160.753606Cr-0.9197260.9021390.698250.6753990.69507Cu-0.983480.502780.669660.621675Hg-0.4791390.6124650.590405Ni-0.577620.634098Pb-0.738756根據(jù)以上數(shù)據(jù)分析，可以得到各重金屬元素大致的空間分布，分別為：1. Hg、Cd、Pb、Zn、Cr均主要分布于工業(yè)區(qū)和左下方的部分交通區(qū)；2. As、Cu主要分布于工業(yè)區(qū)；3. Ni的分布比較分散，沒有明顯的富集現(xiàn)象。（二）重金屬污染程度的劃分為了全面、綜合地反映污染物的整體污染水平

12、，需要一種同時考慮污染物綜合污染水平的多因子評價方法，這里采用內(nèi)梅羅指數(shù)法（內(nèi)梅羅指數(shù)評價法的修正及其應(yīng)用），其計算公式如下：式中，P為土壤污染元素綜合污染指數(shù)；為土壤中污染元素i的實測值；為土壤中污染元素i的評價標(biāo)準(zhǔn)值。內(nèi)梅羅指數(shù)土壤污染評價標(biāo)準(zhǔn)如表二。表二 內(nèi)梅羅指數(shù)土壤污染評價標(biāo)準(zhǔn)等級內(nèi)梅羅污染指數(shù)污染等級IP0.8清潔（安全）II0.7P2.50尚清潔（警戒線）III2.504.25輕度污染IV4.257.20中度污染VP7.20重污染但由于內(nèi)梅羅指數(shù)法過分突出極大值對水質(zhì)污染的影響，會出現(xiàn)評價項目中只有一項指標(biāo)值偏高，而其他指標(biāo)值均較低也會使綜合評分值偏高的情況，故對內(nèi)梅羅法進(jìn)行改進(jìn)

13、如下：其中 式中：為權(quán)重值前n項的平均值，n值根據(jù)評價數(shù)據(jù)確定；為前n項的值；m為前n項中的項數(shù)。對于以上8種重金屬，取權(quán)重值如表三：表三 各重金屬權(quán)重值A(chǔ)sCdCrCuHgNiPbZn權(quán)重值0.0700.1400.0140.0200.6980.0300.0700.020計算得到域內(nèi)不同重金屬的污染程度評價結(jié)果如表四：表四 各區(qū)域內(nèi)重金屬污染程度AsCdCrCuHgNiPbZn生活區(qū)IIIIIIIIIIIIIIIIIIIIV工業(yè)區(qū)IIIIIIIIIVVIIIIIIV山區(qū)IIIIIIIIIIIIIIIII交通區(qū)IIIIIIIVVVIIIIIIIV公園綠地區(qū)IIIIIIIIIIIVIIIIIIV（

14、三）重金屬污染原因分析根據(jù)（一）、（二）中所得數(shù)據(jù)及結(jié)論，將各重金屬污染的原因做如下分析：1. 對于除Ni以外的其他重金屬元素，均有60%以上的重污染點分布在工業(yè)區(qū)，由此可認(rèn)為重金屬污染主要源自工業(yè)區(qū)，這主要是由于工業(yè)區(qū)中的工廠在生產(chǎn)過程中產(chǎn)生了含該種金屬的廢棄物，在排放時也沒有進(jìn)行必要的凈化處理而直接排入水體或大氣中，從而造成了污染的形成與擴(kuò)散，2. Hg、Cd、Pb、Zn四種金屬元素的均有約20%重污染點分布在交通區(qū)，這些污染主要來自于公路上行駛車輛燃燒汽油或柴油所產(chǎn)生的尾氣沉降及輪胎的磨損。3. 金屬元素Zn在生活區(qū)和公園綠地區(qū)也有一定的分布，這可能是由于該區(qū)域所使用的肥料中含有Zn元素

15、。4. 金屬元素Ni的污染點較為分散，可認(rèn)為該污染與人為因素關(guān)系不大，主要來自于土壤母質(zhì)本身。（四）重金屬污染物的傳播特征重金屬一般不會直接通過土壤進(jìn)行傳播，它一般通過如下方式進(jìn)行傳播：1. 空氣傳播 工業(yè)區(qū)中的工廠會排放出大量的工業(yè)廢氣，這些廢氣中可能含有大量的重金屬顆粒，這些顆粒會隨空氣的流動而分散到城市各處，從而造成污染，另外，廢氣也可能會隨降雨滲入土壤。2. 水源傳播 工廠排出的工業(yè)廢水中含有大量的重金屬離子及沉淀物，這些物質(zhì)均可隨水滲入土壤，造成污染。3. 工業(yè)廢渣和生活垃圾 廢渣及生活垃圾的堆積會造成重金屬元素的大量積聚。圖三 城區(qū)地勢圖另外，重金屬的傳播受地形因素的影響較為明顯，

16、從數(shù)據(jù)來看，該城市整體海拔較低，地勢起伏不大。相對1000米的取樣間隔，300米的高度差不會因高空吸附或雨水沖刷造成污染物沉積異常。但山勢對風(fēng)向的影響勢必導(dǎo)致污染物濃度區(qū)域性變化。圖四 Cd污染程度分布圖如上圖的標(biāo)注，山腳、較易形成污染物富集；由于、山谷的存在，經(jīng)過山谷的大風(fēng)在A、B區(qū)域減速后極易產(chǎn)生污染物富集；此外，迎向污染源的山坡受到的污染會明顯大于背向污染源的山坡。從Hg、Cd的污染分布情況可以驗證這幾點（如圖四）：（五）建立模型確定污染源在進(jìn)行污染源的確定之前先進(jìn)行篩選，如果污染是由點污染源引起的，那么污染源就一定處于污染程度最大，即重金屬濃度最高的采樣點的附近。由于城市中建筑物較多，

17、人類活動頻繁，導(dǎo)致局部氣流復(fù)雜，因此通過空氣傳播占主要地位，又由于城市缺少地表水系，雨水及其它污水會通過下水管排出城外，故影響相對較小，因此只討論通過空氣傳播的模型：設(shè)是t時刻點處某種重金屬污染物的濃度。任取一個閉合曲面S，它所圍的區(qū)域是，由于氣體擴(kuò)散，由t到t+時刻的這段時間內(nèi)，通過S流入的質(zhì)量為：由高斯公式得：其中分別是沿方向的擴(kuò)散系數(shù)。由于衰減（例如吸收、代謝等），的質(zhì)量減少為：其中為衰減系數(shù)。由于重金屬污染物在土壤中很難降解，因此衰減系數(shù)可以取為0。由質(zhì)量守恒定律，在內(nèi)由于擴(kuò)散和衰減的合作用，積存于內(nèi)物質(zhì)的總量為從另外的一個角度來看，內(nèi)由于濃度的變化引起的質(zhì)量增加為：顯然，即由設(shè)擴(kuò)散源

18、在點(由假設(shè)可知，如果認(rèn)為經(jīng)過相當(dāng)長的時間后，擴(kuò)散已經(jīng)停止，則物質(zhì)處于平衡狀態(tài)，有：即：由于污染物在空中傳播時，地面上污染物的沉積量與z方向上的傳播無關(guān)，故上述方程可以改寫成：因此問題轉(zhuǎn)化成二維含源反演問題，可以用貝葉斯推理來估計參數(shù)和。貝葉斯推理的基礎(chǔ)是貝葉斯定理,它可以表述如下: 其中,為模型參數(shù)，y為觀測數(shù)據(jù)，為參數(shù)的先驗概率密度數(shù)，為似然函數(shù)，為參數(shù)的后驗概率密度函數(shù)。在貝葉斯推理中，p()表示在未獲得測量數(shù)據(jù)之前，對模型參數(shù)分布的認(rèn)識，主要來源于以往數(shù)據(jù)、 經(jīng)驗和主觀判斷等。代表模型參數(shù)擬合測量數(shù)據(jù)的程度，越大表示擬和效果越好，反之則差。表示了獲得測量數(shù)據(jù)以后模型參數(shù)的分布規(guī)律，即為

19、在統(tǒng)計反演意義下的反問題的解。其步驟通常分為以下三步：1. 確定先驗概率密度數(shù)一般情況下，對于參數(shù)，如果僅知道它分布的區(qū)間，那么可以假設(shè)先驗概率為均勻分布：=假設(shè)參數(shù)之間相互獨立則：其中m為模型參數(shù)的個數(shù)。2. 確定似然函數(shù)似然函數(shù)的構(gòu)造對于貝葉斯推理的結(jié)果具有很大影響，一般可以認(rèn)為測量誤差服從正態(tài)分布，此時似然函數(shù)可以表示為：其中n為測量個數(shù)， 可以用解析解表示，一般情況下，為數(shù)值解。3. 對后驗概率密度函數(shù)抽樣，進(jìn)而獲得參數(shù)的估計值選取馬爾科夫鏈蒙特卡洛法(Markov chain Mont Carlo，簡記為MCMC)。算法為MH算法，表示如下：(1) 設(shè)定參數(shù)初始值并設(shè)定；(2) pr

20、oposal分布中生成，計算接受概率：(3) 生成0,1之間均勻分布的隨機(jī)數(shù)；(4) 如果 ，設(shè)置；否則。當(dāng) (迭代次數(shù))時，停止；否則，繼續(xù)第(2)步。其中，為目標(biāo)概率密度函數(shù)，在貝葉斯推理中為似然函數(shù)。（六）模型的求解通過MATLAB編程，使用MH算法進(jìn)行50000次抽樣，使目標(biāo)坐標(biāo)值趨于平穩(wěn)，最終得到各重金屬元素污染源的地理位置坐標(biāo)如表五。（圖五為Hg元素進(jìn)行抽樣后的迭代曲線圖，可以看到圖線最終趨于一條平穩(wěn)的直線。）圖五 Hg元素抽樣后的迭代曲線圖表五 各重金屬元素污染源地理位置坐標(biāo)AsCdCrCuHgNiPbZnX值27422382347027072708164728824152Y值72943693230922952294272936172299（七）對城市地質(zhì)環(huán)境的演變模式的探究地質(zhì)環(huán)境的演變模式是一個地域在一段較長時期內(nèi)所表現(xiàn)出來的地質(zhì)變化規(guī)律與