全國大學(xué)生數(shù)模競賽論文城市表層土壤重金屬污染分析

1、城市表層土壤重金屬污染分析摘 要隨著城市經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷增加，人類活動(dòng)對城市環(huán)境質(zhì)量的影響日顯突出。本文針對某城市地表土壤重金屬測量數(shù)據(jù)，開展城市環(huán)境質(zhì)量評(píng)價(jià)，研究人類活動(dòng)影響下城市地質(zhì)環(huán)境的演變模式。問題1：對8種重金屬元素進(jìn)行空間散亂點(diǎn)插值和Kriging插值，利用ArcGIS地統(tǒng)計(jì)軟件得到8種重金屬的空間分布圖，運(yùn)用單項(xiàng)污染指標(biāo)法和內(nèi)梅羅綜合污染指標(biāo)法分別得出每一種重金屬在不同功能區(qū)的污染等級(jí)：AsCd CrCuHgNiPbZn一區(qū)輕污染中污染中污染重污染中污染輕污染中污染重污染二區(qū)中污染重污染輕污染重污染重污染輕污染重污染重污染三區(qū)輕污染輕污染輕污染輕污染輕污染輕污染輕污

2、染輕污染四區(qū)輕污染輕污染輕污染重污染重污染輕污染中污染中污染五區(qū)輕污染輕污染輕污染中污染重污染輕污染輕污染輕污染綜合污染等級(jí)如下：功能區(qū)污染指標(biāo)污染等級(jí)生活區(qū)3.4349重污染工業(yè)區(qū)4.0279重污染山區(qū)1.0622輕度污染交通區(qū)3.5196重污染公園綠地區(qū)2.2354中度污染問題2：對8種土壤重金屬元素含量運(yùn)用相關(guān)性分析，主成分分析和聚類分析發(fā)現(xiàn)，該市土壤中重金屬污染主要受鐵路交通，農(nóng)業(yè)農(nóng)藥塑料薄膜使用，污灌，生活垃圾堆放，采礦這六個(gè)方面影響問題3：重金屬的傳播特征一般是離污染源越遠(yuǎn)，污染程度越小。將問題1中空間分布圖對應(yīng)在該區(qū)域的所有插值點(diǎn)作為樣本點(diǎn)，建立層次規(guī)劃模型，目標(biāo)層為整個(gè)區(qū)域的某

3、元素加權(quán)平均濃度，準(zhǔn)則層為該元素污染源，決策層為該元素污染源內(nèi)樣本點(diǎn)的濃度，建立非線性規(guī)劃模型，得到兩組最優(yōu)權(quán)重表示污染源內(nèi)的點(diǎn)在該污染源中所占比例，表示污染源內(nèi)的點(diǎn)對樣本點(diǎn)加權(quán)濃度的影響度。根據(jù)越小，污染源對加權(quán)平均濃度的影響力越大，對As進(jìn)行模擬，得到，可以判斷顯然為污染源，是污染源的可能性很小。越大，其對應(yīng)點(diǎn)的As濃度越大，越接近污染源的中心。問題4：將單一的測定土壤中的重金屬濃度，擴(kuò)展為測定空氣，水域，植物和土層（050cm）的各重要元素濃度，并將取樣間距縮小,定期對其取樣，將m種元素n個(gè)時(shí)期的元素濃度均值作為研究對象，以每次檢測時(shí)各元素污染指數(shù)均值作為參考數(shù)列，進(jìn)行灰色關(guān)聯(lián)度分析，得

4、到各個(gè)時(shí)期主要污染元素的演變；并對每一種元素濃度污染指標(biāo)進(jìn)行灰色預(yù)測，得到各元素的進(jìn)入土壤和在土壤中的傳播規(guī)律及演變模式，結(jié)合地理統(tǒng)計(jì)軟件ArcGIS，畫出各元素在不同介質(zhì)中的空間分布圖，對各圖進(jìn)行比較，分析，全面地分析城市地質(zhì)環(huán)境的演變模式。 關(guān)鍵詞：地統(tǒng)計(jì)學(xué)，內(nèi)梅羅綜合污染評(píng)價(jià)，多元分析，灰色預(yù)測 一、問題重述隨著城市經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷增加，人類活動(dòng)對城市環(huán)境質(zhì)量的影響日顯突出。對城市土壤地質(zhì)環(huán)境異常的查證，以及如何應(yīng)用查證獲得的海量數(shù)據(jù)資料開展城市環(huán)境質(zhì)量評(píng)價(jià)，研究人類活動(dòng)影響下城市地質(zhì)環(huán)境的演變模式，日益成為人們關(guān)注的焦點(diǎn)。按照功能劃分，城區(qū)一般可分為生活區(qū)、工業(yè)區(qū)、山區(qū)、

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

6、求通過數(shù)學(xué)建模來完成以下任務(wù)：(1) 給出8種主要重金屬元素在該城區(qū)的空間分布，并分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。(2) 通過數(shù)據(jù)分析，說明重金屬污染的主要原因。(3) 分析重金屬污染物的傳播特征，由此建立模型，確定污染源的位置。(4) 分析所建立模型的優(yōu)缺點(diǎn)，為更好地研究城市地質(zhì)環(huán)境的演變模式，還應(yīng)收集什么信息？有了這些信息，如何建立模型解決問題？二、模型假設(shè)1. 假設(shè)該城市土壤，地質(zhì)環(huán)境短時(shí)間內(nèi)不發(fā)生變化。2. 假設(shè)采集回的樣本能很好的反應(yīng)采樣點(diǎn)的實(shí)際情況。3. 假設(shè)題目所給數(shù)據(jù)真實(shí)準(zhǔn)確。4. 假設(shè)在土壤采樣，化驗(yàn)，分析過程中未出現(xiàn)任何操作錯(cuò)誤或條件改變。5. 假設(shè)采樣過程中未出現(xiàn)特

7、大嚴(yán)重重金屬污染情況三、問題分析根據(jù)題目所給信息，對城市表層土壤重金屬污染進(jìn)行分析討論后，解決了如下四個(gè)問題： （1）土壤重金屬空間分布是應(yīng)用地統(tǒng)計(jì)學(xué)方法對變量空間結(jié)構(gòu)的表達(dá)與模擬，為了更好的說明土壤中重金屬的分布情況，首先要對數(shù)據(jù)進(jìn)行預(yù)處理，由于在土壤采樣，化驗(yàn)，分析過程中會(huì)出現(xiàn)操作錯(cuò)誤或條件改變，這將會(huì)造成變量Kriging插值連續(xù)表面中斷，使得實(shí)驗(yàn)半方差函數(shù)1發(fā)生畸變，甚至?xí)谏w變量固有的空間結(jié)構(gòu)特征，因此首先要剔除異常值，然后進(jìn)行變量正態(tài)分布檢驗(yàn)，然后計(jì)算半方差函數(shù)獲取土壤變量空間自相關(guān)性，最后用Kriging插值法2獲取土壤變量的空間分布圖；應(yīng)用綜合污染指數(shù)法，能全面綜合地反映土壤的

8、污染程度，所以要分析城區(qū)內(nèi)不同區(qū)域重金屬的污染程度，需采用兼顧單元素污染指數(shù)平均值和最大值的尼梅羅綜合污染指數(shù)法3，計(jì)算公式如下： （2）為了說明重金屬污染的主要原因，首先要對這8種重金屬元素進(jìn)行相關(guān)性分析，以確定它們的污染來源是否相同，或具有關(guān)聯(lián)出現(xiàn)，然后再對該8種重金屬元素進(jìn)行主成分分析，考慮占主導(dǎo)地位的成分的影響，使具有強(qiáng)相關(guān)性的元素歸并為一類。最后在前面的多元分析下說明重金屬污染的主要原因。 （3）重金屬的傳播特征一般是離污染源越遠(yuǎn)，污染程度越小。由于同一區(qū)域土壤中重金屬污染的來源可以是單一的，也可以是多途徑的，因此由第（1）問各元素在不同功能區(qū)濃度的空間分布與第（2）問中對各元素污染

9、原因分析，可以確定各污染源的大概位置以及來自于同一污染源的元素種類(即如果污染點(diǎn)和污染源均相同，則可以得出這些元素來自同一污染源)，得出8種元素的可能污染源位置和來自同一可能污染源的元素種類，在個(gè)可能污染源周圍劃分一個(gè)面積為的區(qū)域，將問題（1）中空間分布圖對應(yīng)在該區(qū)域的所有插值點(diǎn)作為樣本點(diǎn)（以下稱樣本點(diǎn)），建立以整個(gè)區(qū)域的某元素加權(quán)平均濃度為目標(biāo)層，每種元素污染源為準(zhǔn)則層，每種元素污染源內(nèi)樣本點(diǎn)的濃度為決策層的非線性規(guī)劃模型。（4）研究城市地質(zhì)環(huán)境的演變模式是個(gè)復(fù)雜的系統(tǒng)過程，僅憑土壤中的重金屬濃度進(jìn)行分析遠(yuǎn)遠(yuǎn)不夠，為此，我們將單一對土壤中重金屬濃度的測定，擴(kuò)展為測定空氣，水域，植物和土層（0

10、50cm）的各重要元素濃度，并將取樣間距縮小，定期對其取樣，記錄。運(yùn)用地理統(tǒng)計(jì)軟件ArcGIS，畫出各元素在不同介質(zhì)中的空間分布圖，對各圖進(jìn)行比較，分析；將m種元素n個(gè)時(shí)期的元素濃度均值作為研究對象，以每次檢測時(shí)各元素濃度的均值作為參考數(shù)列，進(jìn)行灰色關(guān)聯(lián)度分析；對每一種元素濃度污染指標(biāo)進(jìn)行灰色預(yù)測，分析結(jié)果，得到各元素的進(jìn)入和傳播規(guī)律及這種方式轉(zhuǎn)移情況，并得出各元素土層分布規(guī)律，從而全面地分析了城市地質(zhì)環(huán)境的演變模式。四、符號(hào)說明定義符號(hào)說明半方差函數(shù)樣本間距間距為h的“樣本對”對數(shù)樣本點(diǎn)的空間坐標(biāo)樣本點(diǎn)出的元素濃度值插值點(diǎn)處的估計(jì)值污染物的污染指數(shù)污染物的實(shí)測值的平均值污染物的評(píng)價(jià)標(biāo)準(zhǔn)單項(xiàng)污

11、染指數(shù)法評(píng)價(jià)標(biāo)準(zhǔn)重金屬元素j的污染指數(shù)單元素污染指數(shù)最大值單元素污染指數(shù)平均值普通Kriging系數(shù)樣本數(shù)樣本點(diǎn)的加權(quán)平均濃度元素第i個(gè)污染源t污染源個(gè)數(shù)第i個(gè)污染源中樣本點(diǎn)的個(gè)數(shù)I整個(gè)區(qū)域內(nèi)插值點(diǎn)的個(gè)數(shù)污染源內(nèi)的點(diǎn)對樣本點(diǎn)加權(quán)濃度的影響度污染源內(nèi)的點(diǎn)在該污染源中所占比例名詞解釋：1.城市土壤：指由于人為的、非農(nóng)業(yè)作用形成的，并且由于土地的混合、填埋或污染而形成的厚度大于或等于50cm的城區(qū)或郊區(qū)土壤2.地質(zhì)環(huán)境：指地殼上部包括巖石指地殼上部包括巖石、水、氣和生物在內(nèi)的互相關(guān)聯(lián)的系統(tǒng)五、模型建立與求解5.1問題（1）模型的建立與求解5.1.1重金屬元素空間分布模型的建立我們根據(jù)題意在解決重金屬

12、元素空間分布時(shí)建立了兩個(gè)模型：5.1.1.1模型一：根據(jù)給出數(shù)據(jù)，進(jìn)行散亂點(diǎn)插值，利用matlab做出三維坐標(biāo)圖（程序見附錄一）。一般情況下，x，y，z是維數(shù)相同的矩陣，x，y是網(wǎng)格坐標(biāo)矩陣，z是網(wǎng)格點(diǎn)上的高度矩陣，c用于指定在不同高度下的顏色范圍。c省略時(shí)，Matlab認(rèn)為c=z，也即顏色的設(shè)定是正比于圖形的高度的，但是此題圖像中顏色的深淺用于標(biāo)定金屬元素的濃度。（紅>黃>藍(lán)>綠>紫），圖（1）到圖（8）分別表示As,Cd ,Cr,Cu,Hg,Ni,Pb,Zn在空間的濃度大小。 圖（1） 圖（2） 圖（3） 圖（4） 圖（5） 圖（6） 圖（7） 圖（8）5.1.1.

13、2模型二：（1）為使Kriging插值表面連續(xù)，要剔除異常值，利用matlab（程序見附錄二），首先要對數(shù)據(jù)經(jīng)行預(yù)處理（見附件1）。（2）正態(tài)分布檢驗(yàn) 利用spss數(shù)據(jù)統(tǒng)計(jì)分析功能，通過P-P圖檢驗(yàn)數(shù)據(jù)是否符合指定的分布，當(dāng)數(shù)據(jù)符合指定的分布時(shí)，P-P圖中各點(diǎn)近似呈現(xiàn)一條直線，并且我們對于P-P圖中各點(diǎn)不成直線但有一定規(guī)律的，對其數(shù)據(jù)進(jìn)行轉(zhuǎn)換，使轉(zhuǎn)換后數(shù)據(jù)更接近指定分布。對各重金屬含量概率服從做如下表（1）所列假設(shè)： 表（1）假設(shè)8種重金屬元素服從對數(shù)正態(tài)分布重金屬元素AsCdCrCuHgNiPbZn概率分布對數(shù)對數(shù)對數(shù)對數(shù)對數(shù)對數(shù)對數(shù)對數(shù)得到各重金屬元素P-P圖如下所示：圖（9）8種重金屬元

14、素的P-P圖因此題目所給數(shù)據(jù)符合指定的對數(shù)正態(tài)分布（3）半方差函數(shù)分析首先建立如下所示的半方差函數(shù)： (1)其中是樣本間距，又稱為位差；是間距為的“樣本對”對數(shù)，為空間坐標(biāo)，為處的重金屬濃度。然后我們又利用了Kriging插值，插值函數(shù)如下： (2)其中為處得估計(jì)值，為處得觀測值，為普通Kriging系數(shù)，為觀測點(diǎn)個(gè)數(shù)。模型一是基于三維非網(wǎng)格形式的插值擬合，對于土說明壤中重金屬的空間分布誤差較大，所以我們選擇模型二，Kriging插值，較能說明空間分布。5.1.2重金屬元素空間分布模型的求解結(jié)合多元統(tǒng)計(jì)學(xué)和GIS技術(shù)利用ArcGis地統(tǒng)計(jì)學(xué)軟件對土壤的屬性數(shù)據(jù)進(jìn)行數(shù)學(xué)分析和空間制圖可以很直觀的

15、了解研究城市土壤重金屬元素的來源和空間分布特征如下圖所示（以As元素為例，其它元素分布圖見附錄三）： 圖（10） As元素空間分布圖 5.1.3重金屬污染程度模型的建立與求解污染的評(píng)價(jià)方法：（1） 單項(xiàng)污染指數(shù)法 (3) 其中為污染物的污染指數(shù)，為污染物的實(shí)測值的平均值，即對本題所給附表二中給出的數(shù)據(jù)進(jìn)行處理，為污染物的評(píng)價(jià)標(biāo)準(zhǔn)單項(xiàng)污染指數(shù)法評(píng)價(jià)標(biāo)準(zhǔn)，即本題中附表三給出的背景值：非污染，輕污染，中污染，重污染。（程序見附錄四）表（2）各區(qū)域8種重金屬污染程度 元素 區(qū)域As,Cd ,Cr,Cu,Hg,Ni,Pb,Zn一區(qū)輕污染中污染中污染重污染中污染輕污染中污染重污染二區(qū)中污染重污染輕污染重污

16、染重污染輕污染重污染重污染三區(qū)輕污染輕污染輕污染輕污染輕污染輕污染輕污染輕污染四區(qū)輕污染輕污染輕污染重污染重污染輕污染中污染中污染五區(qū)輕污染輕污染輕污染中污染重污染輕污染輕污染輕污染數(shù)據(jù)分析：將每種金屬元素在五個(gè)區(qū)中污染最嚴(yán)重的用紅色標(biāo)出，可以看出二區(qū)紅色最多，因?yàn)槎^(qū)是工業(yè)區(qū)，這符合實(shí)際。而三區(qū)沒有一個(gè)紅色，因?yàn)槿齾^(qū)是山區(qū)，離污染源最遠(yuǎn)，這也符合實(shí)際。（2）內(nèi)梅羅綜合污染指標(biāo)法單因子污染指數(shù)法只能分別反應(yīng)各個(gè)污染物的污染程度，不能全面，綜合地反應(yīng)土壤的污染程度時(shí)，需將單因子污染指數(shù)按一定方法綜合起來進(jìn)行評(píng)價(jià)，即應(yīng)用綜合污染指數(shù)法評(píng)價(jià)。重金屬元素中和污染評(píng)價(jià)采用兼顧單元污染指數(shù)平均值和最大值的

17、內(nèi)梅羅綜合污染指數(shù)法。計(jì)算公式如下： ， (4) (5)式中： 為監(jiān)測點(diǎn)單項(xiàng)污染指標(biāo)平均值中的最大值，為檢測點(diǎn)單項(xiàng)污染指標(biāo)平均值（是基于單項(xiàng)污染指數(shù)法中的八種元素各自的平均值在對其求和除以八得出的結(jié)果）， 為單項(xiàng)指標(biāo)的背景值，為監(jiān)測點(diǎn)的綜合污染指標(biāo)，為j監(jiān)測點(diǎn)所有污染物單項(xiàng)污染指數(shù)中的最大值與背景值的商；為j檢測監(jiān)測點(diǎn)所有污染物單項(xiàng)污染指數(shù)的平均值與背景值的商。表（3） 土壤內(nèi)梅羅指數(shù)評(píng)價(jià)標(biāo)準(zhǔn) 等級(jí) 內(nèi)梅羅污染指數(shù) 污染等級(jí) 一 清潔（安全） 二 尚清潔（警戒線） 三 輕度污染 四 中度污染 五 重污染根據(jù)公式，計(jì)算五個(gè)區(qū)的的內(nèi)梅羅指數(shù)分別為（程序見附錄五）： 生活區(qū)： 3.4349(重污染)

18、 工業(yè)區(qū)： 4.0279(重污染) 山區(qū)： 1.0622(輕度污染) 交通區(qū)： 3.5196(重污染) 公園綠地區(qū)：2.2354(中度污染) 由數(shù)據(jù)分析可以看出對于污染指標(biāo)工業(yè)區(qū)指標(biāo)值最大，而山區(qū)的污染指標(biāo)值最小，這種情況也符合實(shí)際，各區(qū)各重金屬元素污染的分布圖（見附錄五）5.2問題（2）模型的建立與分析 為說明土壤中重金屬污染的主要原因，先要對數(shù)據(jù)的統(tǒng)計(jì)學(xué)特征進(jìn)行初步了解，通過描述統(tǒng)計(jì)分析和相關(guān)性分析，能進(jìn)一步弄清土壤重金屬含量數(shù)據(jù)的表觀特征，為其污染的主要原因的說明提供一定的依據(jù)和理論支持。各重金屬元素的背景值如下表（4）所示： 表（4）8種重金屬元素的背景值元素平均值標(biāo)準(zhǔn)偏差范圍As (

19、g/g)3.60.91.85.4Cd (ng/g)1303070190Cr (g/g)3191349Cu (g/g)13.23.66.020.4Hg (ng/g)3581951Ni (g/g)12.33.84.719.9Pb (g/g)3161943Zn (g/g)691441975.2.1土壤重金屬含量的描述統(tǒng)計(jì)分析特征描述統(tǒng)計(jì)分析主要研究如何將實(shí)驗(yàn)或調(diào)查得到的大量數(shù)據(jù)簡縮成有代表性的數(shù)據(jù)，使其能客觀，全面地反映這組數(shù)據(jù)的全貌，可以將其所提供的信息充分顯現(xiàn)出來，為下一步統(tǒng)計(jì)分析提供依據(jù)，對該市319個(gè)城市土壤樣品8種重金屬含量進(jìn)行統(tǒng)計(jì)分析，結(jié)果見表（5）表（5）城市土壤重金屬含量統(tǒng)計(jì)描述元素

20、樣本數(shù)最小值最大值均值標(biāo)準(zhǔn)差偏度峰度變異系數(shù)（%）As3191.6130.135.67653.02433.32419.69653.28Cd319401619.8302.3962224.98762.0245.64174.4Cr31915.32920.8453.509770.00189.445104.166130.82Cu3192.292528.4855.0167162.915112.755180.685296.12Hg3198.5716000299.71131629.53988.59674.148540.7Ni3194.27142.517.26189.94147.10981.84457.59Pb

21、31919.68472.4861.740950.05784.40628.30881.08Zn31932.863760.82201.2026339.23256.4953.591168.6 與該市土壤環(huán)境背景值進(jìn)行比較，可以發(fā)現(xiàn)，8種重金屬元素在該地區(qū)的平均含量分別是其自然環(huán)境背景值的1.58,2.33,1.73,4.17,37.46,4.54,10.29,14.07倍，都具有明顯的累積現(xiàn)象；變異系數(shù)（C-V）是表征樣本間變異程度的重要參數(shù)，C-V<10%為弱變異；C-V=10%-30%為中等變異，C-V>30%為強(qiáng)變異。該市城市土壤樣品中8種重金屬含量差異很大（見表（5），變異系數(shù)排

22、序依次為：As<Ni<Cd<Pb<Cr<Zn<Cu<Hg，這些重金屬的變異系數(shù)均高于30%，變異系數(shù)越大說明土壤重金屬污染受人為活動(dòng)干擾越強(qiáng)烈，或者可以理解為污染程度越大，從變異系數(shù)中可以看出，Hg含量變異最為明顯，其次為Cu，Zn。5.2.2土壤重金屬含量之間的統(tǒng)計(jì)分析 研究土壤重金屬含量之間的相關(guān)性，在一定程度上可反映這些元素污染程度的相似性或污染元素有相似的來源，利用spss統(tǒng)計(jì)分析軟件包，對該市城市土壤和重金屬含量做pearson相關(guān)分析，主成分分析及聚類分析。（1）相關(guān)性分析：一般情況下，如果元素元素含量間顯著相關(guān)，說明他們來源相同的可能性較

23、大，這一來源既有可能出自天然，也有可能來自人類活動(dòng)造成的符合污染，如表（6）所示，由表可知，城市土壤重金屬除As-Hg外其余其余都具有較高的相關(guān)性，表明它們之間的同源性很高，而As-Hg的相關(guān)性很弱。表（6）重金屬含量的pearson相關(guān)矩陣（n=319）元素AsCdCrCuHgNiPbZnAs10.255*0.189*0.160*0.0640.317*0.290*0.247*Cd0.255*10.352*0.397*0.265*0.329*0.660*0.431*Cr0.189*0.352*10.532*0.103*0.716*0.383*0.424*Cu0.160*0.397*0.532*

24、10.417*0.495*0.520*0.387*Hg0.0640.265*0.103*0.417*10.103*0.298*0.196*Ni0.317*0.329*0.716*0.495*0.103*10.307*0.436*Pb0.290*0.660*0.383*0.520*0.298*0.307*10.494*Zn0.247*0.431*0.424*0.387*0.196*0.436*0.494*1注：*P0.05水平下，顯著相關(guān)；*P0.01水平上，極顯著相關(guān)（2）主成分分析：對城市土壤中這8種重金屬As，Cd,Cr,Cu,Hg,Ni,Pb和Zn進(jìn)行主成分分析，當(dāng)把因子固定為6個(gè)時(shí)，通

25、過spss進(jìn)行主成分分析時(shí)，這6個(gè)因子所對應(yīng)的累積貢獻(xiàn)率為93.156%（見表（7），因而可以判定這6個(gè)因子都為它的主成分因子，從主成分分析的旋轉(zhuǎn)成分矩陣（見表（8）上可以看出 表（7）6個(gè)因子所對應(yīng)的累積貢獻(xiàn)率解釋的總方差成份初始特征值提取平方和載入旋轉(zhuǎn)平方和載入合計(jì)方差的 %累積 %合計(jì)方差的 %累積 %合計(jì)方差的 %累積 %13.56044.50044.5003.56044.50044.5001.82722.83422.83421.15014.37758.8771.15014.37758.8771.52519.06341.8973.96512.06370.941.96512.06370.

26、9411.04513.06854.9654.7689.59680.537.7689.59680.5371.03212.89867.8635.5787.22087.756.5787.22087.7561.02012.74480.6076.4325.39993.156.4325.39993.1561.00412.54993.156表（8）6個(gè)主成分因子的得分系數(shù)矩陣旋轉(zhuǎn)成份矩陣a成份123456As0.1350.1300.0160.9740.0430.084Cd0.2230.9180.1480.0860.0310.126Cr0.8590.187-0.0160.0040.2450.141Cu0.39

27、50.1770.2730.0310.8100.110Hg0.0170.1370.9670.0160.1730.069Ni0.8900.1010.0510.1960.1030.170Pb0.0730.7180.0610.1620.4960.275Zn0.2610.2350.0820.0970.1220.916提取方法 :主成分分析法。 旋轉(zhuǎn)法 :具有 Kaiser 標(biāo)準(zhǔn)化的正交旋轉(zhuǎn)法。 第一主成分中主要是Ni和Cr，第二主成分中主要是Cd和Pb，第三主成分中Hg占主導(dǎo)作用；第四主成分中As占主導(dǎo)作用；第五主成分中Cu占主導(dǎo)作用，第六主成分中Zn占主導(dǎo)作用。結(jié)合相關(guān)性分析和主成分分析結(jié)果我們可以得

28、出：在主成分因子1中，Ni，Gr的因子負(fù)荷量較高，而且他們之間存在顯著相關(guān)性，這說明它們污染來源可能相同，可能主要來源于工業(yè)污染。在主成分因子2中，Cd和Pb的因子負(fù)荷量較高，可能來源于礦物的開采含鉛汽油的排放。在主成分因子3中，Hg的因子負(fù)荷量較高，占主導(dǎo)作用，除來源于土壤母質(zhì)外，還包括人為因素諸如：工業(yè)生產(chǎn)廢料，城市生活垃圾堆放，農(nóng)業(yè)施肥以及污灌等。在主成分因子4中，As的因子負(fù)荷量最高，占主導(dǎo)作用，除來源于地質(zhì)巖石外，還有農(nóng)業(yè)化肥，礦山開采，礦石冶煉等人為因素影響，且可能由于點(diǎn)污染。在主成分因子5中，Cu的因子負(fù)荷量較高，占主導(dǎo)作用，可能由于銅礦開發(fā)產(chǎn)生的尾石，礦石，三廢的排放，有機(jī)肥，

29、含銅農(nóng)業(yè)化學(xué)物質(zhì)等，含鉛汽油的排放。在主成分因子6中，Zn占主導(dǎo)作用，與其它幾種重金屬含量關(guān)系較小可能與其獨(dú)特的積累特征如點(diǎn)源污染等有關(guān)。綜上所述，可以推斷出該市城市土壤重金屬污染與工業(yè)和城市生活垃圾污染物的積累，農(nóng)業(yè)，城市交通有關(guān)。（三）聚類分析：然后采用R型聚類分析對八種元素進(jìn)行相關(guān)性處理：使用的方法是最短距離法 它的直觀意義為兩個(gè)類中最近兩點(diǎn)間的距離，R型聚類樹形圖圖見下： 圖（11）重金屬元素聚類樹形圖注：樹狀圖的的橫坐標(biāo)1，2, 3，4,5,6,7，8分別代表As,Cd ,Cr,Cu,Hg,Ni,Pb,Zn這8種元素。得出的結(jié)論與主成分分析類似，因此我們這樣總結(jié)：該市城市土壤重金屬污

30、染與工業(yè)和城市生活垃圾污染物的積累，農(nóng)業(yè)，城市交通有關(guān)。5.3問題（3）模型的建立與求解我們建立了兩個(gè)模型5.3.1模型一的建立與求解由分析根據(jù)問題（1），（2）的結(jié)果可初步判斷：Cd有三個(gè)污染源，Pb有兩個(gè)污染源，其中兩個(gè)污染源是共同的；Cr和Ni是同一污染源；As有兩個(gè)污染源；Cu有一個(gè)污染源；Hg有三個(gè)污染源；Zn有兩個(gè)污染源。設(shè)為某元素第i個(gè)污染源，設(shè)為第i個(gè)污染源的面積，有第一問可知內(nèi)插值點(diǎn)的個(gè)數(shù)，內(nèi)第j點(diǎn)的坐標(biāo)記為，元素濃度記為，以內(nèi)的插值點(diǎn)作為樣本點(diǎn)，對于同一污染源的N種元素其內(nèi)樣本點(diǎn)的坐標(biāo)取值相同，繪制如下層次分析流程圖：目標(biāo)層準(zhǔn)則層決策層 加權(quán)平均濃度污染源污染源污染源樣本點(diǎn)

31、樣本點(diǎn)點(diǎn)樣本點(diǎn)樣本點(diǎn)樣本點(diǎn)樣本點(diǎn) 因此我們建立了非線性規(guī)劃模型：目標(biāo)函數(shù)： (6)約束條件： (7)又上可得到最優(yōu)權(quán)值和，為可能污染源對加權(quán)平均濃度的影響力，越小，污染源對加權(quán)平均濃度的影響力越大，我們可以設(shè)定一個(gè)閥值，若<0.1，則可判斷是污染源；表示第i個(gè)可能污染源中第j個(gè)樣本點(diǎn)在該污染源中所占的比重，顯然，越大，該點(diǎn)越接近于污染源中心。利用Lingo軟件對模型進(jìn)行求解（程序見附錄六）：對As進(jìn)行模擬，得到，可以判斷顯然為污染源，是污染源的可能性很小，區(qū)域濃度差異可能是As的沉積引起的。觀察所有，可知越大，其對應(yīng)點(diǎn)的As濃度越大，也可理解為越接近傳染源的中心，顯然，這是與As的空間分

32、布相符的。同理可得到其它各重金屬元素的污染源的位置。5.3.2模型二的建立與求解（1）.區(qū)域網(wǎng)格搜索法：假設(shè)在指定的區(qū)域內(nèi)只由一個(gè)污染源頭造成的，然后對于指定的區(qū)域內(nèi)在坐標(biāo)下求出金屬濃度的最大值，作為最大污染指標(biāo)（），然后根據(jù)坐標(biāo)做成合適的方陣，將此方陣作為指標(biāo)方陣。接著將城區(qū)插值形成若干個(gè)與指標(biāo)方陣大小相同的方陣，把金屬濃度作為Z軸。將指標(biāo)方陣中的污染指數(shù)與城區(qū)的污染指數(shù)做方差求和，式子如下：在對城區(qū)中的每個(gè)方陣區(qū)域的進(jìn)行比較，找出最小的即為金屬濃度與指標(biāo)最為接近的。因?yàn)橹笜?biāo)是污染最為嚴(yán)重的，所以最小的的區(qū)域認(rèn)為是污染源頭。 第二問中通過聚類分析已經(jīng)將是Ni和Cr分成一類， Cd和Pb分成一

33、類， Hg，As，Cu，Zn分別為一類。基于MATLAB求出金屬經(jīng)過分類的和所對應(yīng)的區(qū)域：金屬Ni和CrCd和PbHg AsCuZn1.90E+063.51E+052.17E+098.05E+031.67E+083.12E+08對應(yīng)區(qū)域(上下代表X，Y軸)（3000,4000）（5000,6000）（11000,12000）（4000,5000）（2000,3000）（2000,3000）（3000,4000）（8000,9000）（3000,4000）（2000,3000）（2000,3000）（2000,3000）拿出比較有代表的As進(jìn)行研究，根據(jù)第一問已經(jīng)得出有金屬濃度的因素決定顏色的深

34、淺，在第一問中的三維圖進(jìn)行變換視角，和As得出數(shù)據(jù)進(jìn)行驗(yàn)證：從圖像中可以看出在此處的得出的數(shù)據(jù)符合第一問中給出的值，所以模型正確。5.4問題（4）模型的分析，建立與求解5.4.1 首先對所建立的模型進(jìn)行分析。模型的優(yōu)點(diǎn)：運(yùn)用多種插值方法進(jìn)行比較，選擇較優(yōu)的插值方法，減小了由于方法不當(dāng)引起的系統(tǒng)誤差。能比較直觀地分析土壤，重金屬的濃度變化。運(yùn)用Kriging插值法，誤差小，與實(shí)際吻合度高，且外推能力強(qiáng)。比較同一重金屬元素在不同區(qū)域的濃度分布，能直觀地反映某一重金屬元素在不同區(qū)域的污染程度。用內(nèi)梅羅綜合污染指標(biāo)法綜合地反應(yīng)土壤的污染程度，解決了單因子污染指數(shù)法的不全面性。模型的缺點(diǎn)：相鄰采樣點(diǎn)距離

35、較大，運(yùn)用Kriging插值可能存在誤差。運(yùn)用Kriging插值時(shí)未考慮海拔影響。只考慮了重金屬在土壤中的傳播，忽略了河流，大氣，其它元素對地質(zhì)環(huán)境的影響。5.4.2需要收集的信息：將城區(qū)劃分為間距為1公里左右的網(wǎng)格子區(qū)域，按照每平方公里隨機(jī)均勻的采取5個(gè)樣本點(diǎn)，再對每個(gè)樣本點(diǎn)的表層土（0-50cm每隔10cm取一個(gè)樣品）進(jìn)行取樣，編號(hào)，并用GPS記錄采樣點(diǎn)位置，記錄所含元素濃度數(shù)據(jù)。對該城市各城區(qū)主要河流進(jìn)行采樣，編號(hào)，并記錄位置和元素含量數(shù)據(jù)。對該城市各個(gè)功能區(qū)的空氣進(jìn)行各元素的濃度測量，記錄采樣點(diǎn)位置，元素含量值，并編號(hào)。對該城市各區(qū)的植物進(jìn)行采樣，編號(hào)，并用GPS記錄采樣點(diǎn)位置，并記錄

36、植物中元素含量值。由于城市土壤的巖石中所含元素是經(jīng)年累月的結(jié)果，需考慮在采樣期間是否進(jìn)行過大面積的土地翻新對同一采樣點(diǎn)利用同一方法定期對上面5個(gè)采樣種類進(jìn)行采樣，記錄濃度，時(shí)間。5.4.3模型的建立：1.運(yùn)用Kriging插值法分別作出每一種元素在地表，水域，大氣的空間分布圖，由圖像得出各元素在地表，水域，大氣中的濃度分布情況，得出各統(tǒng)計(jì)量。2.對該城市城區(qū)采集的植物中各元素濃度進(jìn)行Kriging插值，并繪制其空間分布圖3.對該市空氣中元素濃度進(jìn)行Kriging插值，并繪制空間分布圖，得出統(tǒng)計(jì)量。4.繪制個(gè)各元素在土壤中的垂直方向上的變化圖，分析各元素在各土層的濃度變化。5.每次定期取樣均進(jìn)行

37、以上4個(gè)過程。以定期取樣所得統(tǒng)計(jì)量的均值作為原始數(shù)據(jù)，假設(shè)做了N次定期檢測，第j種元素第i次的濃度均值為。（1） 關(guān)聯(lián)度分析設(shè)影響地質(zhì)環(huán)境演變的元素有m種，將每次檢測的各元素濃度均值與背景值P（其余元素的背景值的確定方法與題給雷同）的比值作為參考數(shù)列其中k表示時(shí)刻。假設(shè)有m個(gè)比較數(shù)列則比較數(shù)列對參考數(shù)列在k時(shí)刻的關(guān)聯(lián)系數(shù)為： （12）其中為分辨系數(shù)比較數(shù)列對參考數(shù)列關(guān)聯(lián)度為： （13）由上可得出m種元素的關(guān)聯(lián)度，關(guān)聯(lián)度越大，說明在污染物中所占的比重就越大。（二）灰色預(yù)測（1） 數(shù)據(jù)檢驗(yàn)與處理對每一種元素進(jìn)行灰色預(yù)測，設(shè)參考數(shù)據(jù)為（表示第k次檢測土壤中某元素濃度的平均值與背景值的比值），計(jì)算數(shù)列

38、的級(jí)比： (14)如果所有的級(jí)比都落在可容覆蓋內(nèi)，則數(shù)列可以作為模型GM（1,1）的數(shù)據(jù)進(jìn)行灰色預(yù)測。否則，需要對數(shù)列作必要的變換處理，使其落入可容覆蓋內(nèi)。即取適當(dāng)?shù)某?shù)c，作平移變換： (15)則使數(shù)列的級(jí)比（2） 建立模型建立模型GM（1,1），對已知參考數(shù)列做一次累加生成數(shù)列： (16)記取的加權(quán)均值，則， 為確定的參數(shù)，記 的白化微分方程模型為： (17)其中a是發(fā)展灰度，b是內(nèi)生控制灰度。由于，取為灰導(dǎo)數(shù)，為背景值，則方程（17）相應(yīng)的灰微分方程為或 (18)即矩陣形式為，其中用最小二乘法球的參數(shù)估計(jì)值為 于是方程（18）有響應(yīng)（特解）4則可以得到預(yù)測值 (19)（3） 檢驗(yàn)預(yù)測值殘

39、差檢驗(yàn)：令殘差為，計(jì)算 (20)如果<0.2,則可認(rèn)為達(dá)到一般要求；如果<0.1,則認(rèn)為達(dá)到較高要求。級(jí)比偏差值檢驗(yàn)首先由參考數(shù)據(jù)，計(jì)算出級(jí)比，再用發(fā)展系數(shù)求出相應(yīng)的級(jí)比偏差 (21)如果<0.2,則可認(rèn)為達(dá)到一般要求；如果<0.1，則認(rèn)為達(dá)到較高要求。（3） 預(yù)測預(yù)報(bào)由GM（1,1）所得到的指定時(shí)區(qū)內(nèi)的預(yù)測值，實(shí)際問題的需要，給出相應(yīng)的預(yù)測預(yù)報(bào)。參考文獻(xiàn)1秦耀東，土壤空間變異研究中的半方差問題，2011/9/10234孫洪泉，地質(zhì)統(tǒng)計(jì)學(xué)及其應(yīng)用，徐州；中國礦業(yè)大學(xué)出版社，1990.91935雷志棟，楊詩秀，田間土壤水分滲透的空間分布，水力學(xué)報(bào)，1987（3）；196劉

40、會(huì)燈，朱飛，MATLAB編程基礎(chǔ)與典型應(yīng)用，北京：人民郵電出版社，2008.77姜啟源，謝金星，葉俊，數(shù)學(xué)模型，北京：高等教育出版社，2003.8附錄附錄一：load F:/2.txtload F:/1.txtnum=X2;x1=num(:,1);x2=num(:,2);x3=num(:,3);x4=num(:,4);x5=num(:,5);x6=num(:,6);x7=num(:,7);x8=num(:,8);y1=X1(:,1);y2=X1(:,2);y3=X1(:,3);y4=X1(:,4);z1=0:300:29000;z2=0:300:19000;z1,z2=meshgrid(z1,

41、z2);z3=griddata(y1,y2,y3,z1,z2,'v4');z4=griddata(y1,y2,x1,z1,z2,'v4');mesh(z1,z2,z3,z4);z4=griddata(y1,y2,x2,z1,z2,'v4');mesh(z1,z2,z3,z4); z4=griddata(y1,y2,x3,z1,z2,'v4');mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x4,z1,z2,'v4');mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x5

42、,z1,z2,'v4');mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x6,z1,z2,'v4');mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x7,z1,z2,'v4');mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x8,z1,z2,'v4');mesh(z1,z2,z3,z4)1.7417932.2304722.22643.7426652.6583051.4912422.2292383.4349082.0142753.0239321.7228769.

43、66180618.353011.6107053.0013174.0279351.1233591.171691.2567641.3119151.1701731.2564061.1792231.0622351.5855682.7693421.8727074.71325212.766361.4322852.049493.5196361.7399212.1580221.4076132.2872513.2854781.2430661.9583412.235395附錄二：x1=5.6765 0.3023 53.509 55.106 0.29971 17.2618 61.740 201.20;x2=3 0.

44、2246 69.9 162.7 1.627 9.9 50 338.7;x3=14.67650.9761263.2090543.2060 5.00046.9618211.74001217.3000x4=3.3235 0.3715 156.1910 432.9940 4.5813 12.4382 88.2600 814.9000load F:/2.txtnum=X2;for i=1:319if num(i,1)<x3(1)&&num(i,1)>x4(1)num(i,1)=num(i,1);else if num(i,1)>2*x1(1) num(i,1)=0;el

45、se if num(i,1)>(1/3)*x1(1) num(i,1)=0; else if num(i,1)<x4(1)&&num(i,1)> (1/3)*x1(1) num(i,1)=x1(1)-3*x3(1); else if num(i,1)>x3(1)&&num(i,1)<2*x1(1) num(i,1)=x1(1)+3*x3(1); end end end endendend附錄三： Cd元素空間分布圖 Cr元素空間分布圖 Cu元素空間分布圖 Hg元素空間分布圖 Ni元素空間分布圖 Pb元素空間分布圖 Zn元素空間分布圖

46、附錄四：num=xlsread('F:/1.xls');for i=1:319num(i,6)=num(i,6)./1000;num(i,9)=num(i,9)./1000;endx=num(1:44,5:12);x1=num(45:80,5:12);x2=num(81:146,5:12);x3=num(147:284,5:12);x4=num(285:319,5:12);x5=3.6000 0.130 31.0000 13.2000 0.035 12.3000 31.0000 69.0000;for i=1:8t1(i)=mean(x(:,i)/x5(i);t2(i)=mea

47、n(x1(:,i)/x5(i);t3(i)=mean(x2(:,i)/x5(i);t4(i)=mean(x3(:,i)/x5(i);t5(i)=mean(x4(:,i)/x5(i);endw1=sqrt(max(t1(i)2+mean(t1(i)2)/2);w2=sqrt(max(t2(i)2+mean(t2(i)2)/2);w3=sqrt(max(t3(i)2+mean(t3(i)2)/2);w4=sqrt(max(t4(i)2+mean(t4(i)2)/2);w5=sqrt(max(t5(i)2+mean(t5(i)2)/2);附錄五： As元素在各區(qū)的分布 Cd元素在各區(qū)的分布 Cr元素

48、在各區(qū)的分布 Cu元素在各區(qū)的分布 Hg元素在各區(qū)的分布 Ni元素在各區(qū)的分布 Pb元素在各區(qū)的分布 Zn元素在各區(qū)的分布附錄六：model:sets:x/1.10/;y/1.10/;com(x,y):w1,w2,z1,z2,z3;endsetsmin=sum(com(i,j):(z0-z3(i,j)2);z0=w11*sum(com(i,j):w1(i,j)*z1(i,j)+w22*sum(com(i,j):w2(i,j)*z2(i,j);w11+w22=1;sum(com(i,j):w1(i,j)=1;sum(com(i,j):w2(i,j)=1;data:z1=0.2269920.6174920.9915621.3441321.6693621.9603922.2089922.4052422.5370522.5897220.4393320.8006221.1417421.4570721.7400821.9831622.1773322.311922.3742122.3492920.5938620.9232321.2287821.5044521.7432621.9371122.0765722.1506822.1466822.049920.6948520.9905421.2590321.493982