城表層土壤重金屬污染分析

2011高教社杯全國大學(xué)生數(shù)學(xué)建模競賽承諾書我們仔細閱讀了中國大學(xué)生數(shù)學(xué)建模競賽的競賽規(guī)則.我們完全明白，在競賽開始后參賽隊員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問題。我們知道，抄襲別人的成果是違反競賽規(guī)則的，如果引用別人的成果或其他公開的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻的表述方式在正文引用處和參考文獻中明確列出。我們鄭重承諾，嚴(yán)格遵守競賽規(guī)則，以保證競賽的公正、公平性。如有違反競賽規(guī)則的行為，我們將受到嚴(yán)肅處理。我們參賽選擇的題號是（從A/B/C/D中選擇一項填寫）：A我們的參賽報名號為（如果賽區(qū)設(shè)置報名號的話）：所屬學(xué)校（請?zhí)顚懲暾娜和瑵髮W(xué)參賽隊員（打印并簽名）：1.趙云波劉雄飛指導(dǎo)教師或指導(dǎo)教師組負責(zé)人（打印并簽名）：日期：2011年9月11日賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）：

2011高教社杯全國大學(xué)生數(shù)學(xué)建模競賽編號專用頁賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）：賽區(qū)評閱記錄（可供賽區(qū)評閱時使用）：評閱人評分備注全國統(tǒng)一編號（由賽區(qū)組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進行編號）城市表層土壤重金屬污染分析摘要：本文中，城市表層土壤金屬污染分析需要綜合不同區(qū)域，不同金屬的綜合影響，根據(jù)隨機的數(shù)據(jù)采樣點，通過統(tǒng)計與插值的分析方法進行處理。首先可以考慮對采樣點進行網(wǎng)格化的數(shù)據(jù)處理，然后通過Kriging方法進行空間散亂點的插值處理。通過對函數(shù)的插值結(jié)果的觀察，與原始采樣數(shù)據(jù)的空間分布相比較，可以發(fā)現(xiàn)二者具有較好的吻合度，結(jié)果令人滿意。對城區(qū)內(nèi)不同區(qū)域重金屬污染程度的綜合評價，要對各重金屬評價指標(biāo)分別加權(quán)。利用熵權(quán)法來確定各重金屬評價指標(biāo)的權(quán)重系數(shù)，熵權(quán)法的實際意義在這里體現(xiàn)得尤為明顯，根據(jù)熵權(quán)法得到的相關(guān)系數(shù)均為正值，這一點也驗證了熵權(quán)法在尋找個金屬污染物權(quán)重時的正確性，然后由綜合權(quán)重進行線性加和，得到各個區(qū)域的綜合評定指標(biāo)，同時根據(jù)金屬含量背景值進行等級標(biāo)準(zhǔn)的劃分，從而確定不同區(qū)域的污染程度，結(jié)果與實際完全符合，說明熵權(quán)法的運用是正確的，從而找到重金屬污染的主要原因。區(qū)號區(qū)類型1生活區(qū)2工業(yè)區(qū)3山區(qū)4主干道5公園綠地污染等級輕度污染重度污染無污染重度污染輕度污染在建立重金屬污染的傳播特征模型，先假設(shè)了污染源的位置，然后考慮根據(jù)擴散定律建立模型，根據(jù)一維擴散方程建立模型，但是這樣的話在數(shù)據(jù)處理上必須首先根據(jù)采樣點濃度特征大致確定污染源的位置，然后建立方程，根據(jù)采樣數(shù)據(jù)提取信息求解，由于采樣點較多，本問題的處理可能會遺漏部分有用信息。如果與實際偏差較大，沒有得到有效結(jié)果，可以進一步尋求其他解決方法。對于問題四，可以基于問題一二三得到的模型進行綜合考慮，從數(shù)據(jù)的處理方式，以及地質(zhì)演變過程中時間變量的影響等因素入手，尋求更好的建立模型的方式。若引入時間變量，將擴散模型修改為：F3,J)+k3,j)(*+*)—f3,j,u)=*ox2ay2dt如此一來，還需要知道濃度的時間變化率，即不同時刻采樣點的濃度數(shù)值。然后要考慮二維（或三維）變差函數(shù)。但它們都是以一維函數(shù)為基礎(chǔ)的，只不過要考慮各向同性或各向異性，以及結(jié)構(gòu)的套合?；蛘呖梢赃M一步的加強考慮海拔Z的影響，得到四元的方程模型，然后建立求解。TOC\o"1-5"\h\z\o"CurrentDocument"一問題重述3\o"CurrentDocument"二、問題分析4\o"CurrentDocument"三模型假設(shè)5\o"CurrentDocument"四變量與符號說明5\o"CurrentDocument"五模型建立與求解6\o"CurrentDocument"1隨機場和區(qū)域化變量6\o"CurrentDocument"2變差函數(shù)的構(gòu)造6\o"CurrentDocument"3平穩(wěn)性假設(shè)和本征假設(shè)6\o"CurrentDocument"4實驗變差函數(shù)7\o"CurrentDocument"5變差函數(shù)的理論模型.76Kriging方法插值8\o"CurrentDocument"7確定各重金屬評價指標(biāo)的權(quán)重系數(shù)138傳播特征16\o"CurrentDocument"六、模型評價與改進18\o"CurrentDocument"1模型評價18\o"CurrentDocument"2改進方向18\o"CurrentDocument"參考文獻19\o"CurrentDocument"附件20\o"CurrentDocument"附件一重金屬空間分布模型的建立與求解的matlab實現(xiàn)20\o"CurrentDocument"附件二利用熵權(quán)法進行各區(qū)域污染程度的綜合評定的matlab實現(xiàn)22一問題重述隨著城市經(jīng)濟的快速發(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ū)域環(huán)境受人類活動影響的程度不同?，F(xiàn)對某城市城區(qū)土壤地質(zhì)環(huán)境進行調(diào)查。為此，將所考察的城區(qū)劃分為間距1公里左右的網(wǎng)格子區(qū)域，按照每平方公里1個采樣點對表層土(0~10厘米深度)進行取樣、編號，并用GPS記錄采樣點的位置。應(yīng)用專門儀器測試分析，獲得了每個樣本所含的多種化學(xué)元素的濃度數(shù)據(jù)。另一方面，按照2公里的間距在那些遠離人群及工業(yè)活動的自然區(qū)取樣，將其作為該城區(qū)表層土壤中元素的背景值。附件1列出了采樣點的位置、海拔高度及其所屬功能區(qū)等信息，附件2列出了8種主要重金屬元素在采樣點處的濃度，附件3列出了8種主要重金屬元素的背景值?，F(xiàn)要求通過數(shù)學(xué)建模來完成以下任務(wù)：給出8種主要重金屬元素在該城區(qū)的空間分布，并分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。通過數(shù)據(jù)分析，說明重金屬污染的主要原因。分析重金屬污染物的傳播特征，由此建立模型，確定污染源的位置。分析你所建立模型的優(yōu)缺點，為更好地研究城市地質(zhì)環(huán)境的演變模式，還應(yīng)收集什么信息？有了這些信息，如何建立模型解決問題？二、問題分析本題中，城市表層土壤金屬污染分析需要綜合不同區(qū)域，不同金屬的綜合影響，根據(jù)隨機的數(shù)據(jù)采樣點，考慮通過統(tǒng)計與插值的分析方法進行處理。對于問題一，需要給出8種主要重金屬的空間分布，并分析城區(qū)不同區(qū)域的污染程度。由于采樣數(shù)據(jù)給出的采樣點相對來說空間位置比較散亂，首先可以考慮對此進行網(wǎng)格化的數(shù)據(jù)處理，然后通過Kriging方法進行空間散亂點的插值處理。Kriging是一種距離加權(quán)的插值方法，從地質(zhì)統(tǒng)計學(xué)中借鑒而來，在空間數(shù)據(jù)場滿足給定的統(tǒng)計分布特征的前提假設(shè)下進行插值。通過對函數(shù)的插值結(jié)果的觀察，與原始采樣數(shù)據(jù)的空間分布相比較，進一步的確定該方法可以獲得的插值效果是否有效。對城區(qū)內(nèi)不同區(qū)域重金屬污染程度的綜合評價，要對各重金屬評價指標(biāo)分別加權(quán)。權(quán)重是衡量因子集中某一因子對土壤污染程度影響相對大小的量,權(quán)重系數(shù)越大，則該因子對土壤的影響程度越大，考慮利用熵權(quán)法來確定各重金屬評價指標(biāo)的權(quán)重系數(shù)，在信息論中信息熵表示系統(tǒng)的有序程度，一個系統(tǒng)的有序程度越高，則信息熵越大，反之，一個系統(tǒng)的無序程度越高，則信息熵越小。所以，可以根據(jù)各項指標(biāo)的指標(biāo)值的差異程度，利用信息熵的這個工具計算出各指標(biāo)的權(quán)重。然后由綜合權(quán)重進行線性加和，可以得到各個區(qū)域的綜合評定指標(biāo)，同時根據(jù)金屬含量背景值進行等級標(biāo)準(zhǔn)的劃分，從而確定不同區(qū)域的污染程度。對于問題二，可以結(jié)合問題一中得到的模型，同時對采樣數(shù)據(jù)進行簡單的分析，根據(jù)各區(qū)域的污染程度的不同，可以感性的得到重金屬污染的主要原因必定和污染最嚴(yán)重的區(qū)域有直接關(guān)系。對于問題三，需要建立重金屬污染的傳播特征模型，并確定污染源的位置。可以考慮根據(jù)擴散定律建立模型，首先借鑒一維擴散方程建立模型，但是這樣的話在數(shù)據(jù)處理上必須首先根據(jù)采樣點濃度特征大致確定污染源的位置，然后建立方程，根據(jù)采樣數(shù)據(jù)提取信息求解，由于采樣點較多，本問題的處理可能會遺漏部分有用信息。如果與實際偏差較大，可以進一步尋求其他解決方法。對于問題四，可以基于問題一二三得到的模型進行綜合考慮，從數(shù)據(jù)的處理方式，以及地質(zhì)演變過程中時間變量的影響等因素入手，尋求更好的建立模型的方式。三模型假設(shè)金屬污染濃度場按穩(wěn)定場處理，即各坐標(biāo)點濃度不隨時間變化。假設(shè)各區(qū)域的土壤特性相同。重力對金屬污染的影響忽略。四變量與符號說明變量符號符號說明Xi取樣點橫坐標(biāo)yi取樣點縱坐標(biāo)『街h)點x處半差函數(shù)y*O)h］的觀測值也各點對應(yīng)權(quán)重系數(shù)q.o估計值點與i點之間的變差函數(shù)值i點與j點之間的變差函數(shù)值某金屬濃度值矩陣某金屬濃度值歸一化矩陣第i個因素下第j個評價值的比重ei第i個因素的熵值W各金屬權(quán)重矩陣d.j第j區(qū)的加權(quán)綜合濃度指標(biāo)五模型建立與求解模型需要給出8種主要重金屬的空間分布，并分析城區(qū)不同區(qū)域的污染程度。由于采樣數(shù)據(jù)給出的采樣點相對來說空間位置比較散亂，首先可以考慮對此進行網(wǎng)格化的數(shù)據(jù)處理，然后通過Kriging方法進行空間散亂點的插值處理。Kriging是一種距離加權(quán)的插值方法，從地質(zhì)統(tǒng)計學(xué)中借鑒而來，在空間數(shù)據(jù)場滿足給定的統(tǒng)計分布特征的前提假設(shè)下進行插值。利用Kriging方法對各濃度散亂點進行插值處理，建立起濃度分布的空間曲面(曲面的高度值即用來表征重金屬濃度值)，具體介紹如下：Kriging方法就是對空間數(shù)據(jù)進行加權(quán)插值的權(quán)值設(shè)計方法。Kriging方法通過引進以距離為自變量的變差函數(shù)來計算權(quán)值。由于變差函數(shù)既可以反映變量的空間結(jié)構(gòu)特性，又可以反應(yīng)變量的隨機分布特性，所以利用Kriging方法進行空間數(shù)據(jù)插值往往可以取得理想的效果。另外，通過設(shè)計變差函數(shù)，Kriging方法很容易實現(xiàn)局部加權(quán)插值，這樣就克服了一般距離加權(quán)插值方法插值結(jié)果的不穩(wěn)定性。1隨機場和區(qū)域化變量首先，重金屬的濃度數(shù)據(jù)場可表示為分布于空間的單值函數(shù)S=f(x，y，z)，由題意知S為標(biāo)量，則數(shù)據(jù)場為標(biāo)量場。運用統(tǒng)計學(xué)的方法來研究該數(shù)據(jù)場，首先將f看成隨機函數(shù)，記為Z，依賴于多個自變量的隨機函數(shù)，稱為隨機場。以空間點x的直角坐標(biāo)為自變量的隨機場稱為一個區(qū)域化變量。區(qū)域化變量在觀測前，可以看作是隨機場；觀測后就得到隨機場的一個實現(xiàn)(一般都記作Z(x)，寫法上不加區(qū)別)。濃度區(qū)域化變量同時反映地質(zhì)變量的結(jié)構(gòu)性和隨機性特征。從地質(zhì)學(xué)的觀點來看，區(qū)域化變量可反映地質(zhì)變量的以下特征：局部性、連續(xù)性、異向性、可遷性。2變差函數(shù)的構(gòu)造假設(shè)空間點x只在一維x軸上變化，我們把區(qū)域化變量Z(x)在x和x+h(h為與x具有相同維數(shù)的距離向量)2個點處的值之差的方差之半定義為Z(x)在x方向上的變差函數(shù)，記為y(x，h)，艮口y(x,h)=-|Var[Z(x)—Z(x-|-]o進一步的，由于點x和h是在二維(或三維)空間中變化的，所以要考慮二維(或三維)變差函數(shù)。但它們都是以一維變差函數(shù)為基礎(chǔ)的，只不過要考慮各向同性或各向異性，還要考慮結(jié)構(gòu)的套合。這里暫時先以各向同性進行數(shù)據(jù)分析處理。3平穩(wěn)性假設(shè)和本征假設(shè)當(dāng)區(qū)域化變量Z(x)滿足下列條件時，則稱Z(x)滿足二階平穩(wěn)(或弱平穩(wěn))。在整個研究區(qū)域內(nèi)，區(qū)域化變量Z(x)的數(shù)學(xué)期望存在，且等于常數(shù)，即E[Z(x)]=m(常數(shù))，Px；在整個研究區(qū)域內(nèi)，區(qū)域化變量Z(x)的協(xié)方差函數(shù)存在且相同(即只依賴于滯后*，而與X無關(guān))，即Cov{Z(k),"x+h)}=E[Z(x)-Z(k+h)]—m2=C(h),Vx,Vho在實際工作中經(jīng)常連二階平穩(wěn)假設(shè)也不能滿足，故提出本征假設(shè)。當(dāng)區(qū)域化變量Z(x)的增量[Z(x)-Z(x+h)]滿足下列條件時，則稱Z(x)滿足本征假設(shè)，或說Z(x)是本征的。在整個研究區(qū)域內(nèi)有，區(qū)域化變量Z(x)的數(shù)學(xué)期望存在，且等于常數(shù)，即E[Z(x)-Z(x+h)]=0，PX，Ph；在整個研究區(qū)域內(nèi)，增量[Z(x)-Z(x+h)]的方差函數(shù)存在且平穩(wěn)(不依賴于X)，即Var[Z(x)—Z(X+h)]=E[Z(x)—Z(x+h)2y(hlVx.Vh4實驗變差函數(shù)實驗變差函數(shù)就是根據(jù)觀測數(shù)據(jù)構(gòu)造變差函數(shù)y(h)的估計值少】。有了二階平穩(wěn)假設(shè)或本征假設(shè)，重金屬濃度區(qū)域化變量Z(x)的增量[Z(x)-Z(x+h)]只依賴于分隔它們的h，而不依賴于具體位置X。這樣，被向量h分割的每一個數(shù)據(jù)對{Z(xi),Z(xi+h)}(i=1，2，，N(h))可以看成是{Z(x),Z(x+h)}一次不同的實現(xiàn)(N(h)是被向量h相隔的數(shù)據(jù)對的個數(shù))。這樣便可以用求所有這些觀測值的算術(shù)平均的方法來計算:：門E于是得到這就是實驗變差函數(shù)的基本計算公式。5變差函數(shù)的理論模型為了對區(qū)域化變量的未知值作出估計，需要將實驗變差函數(shù)擬合成相應(yīng)的理論變差函模型，這些模型將直接參與Kriging方法計算。變差函數(shù)模型可以分為有基臺值和無基臺值兩大類，這里利用常用的有基臺值模型進行計算。3種常用的有基臺值模型如下：(1)球狀模型(亦稱馬特隆模型，在原點處為線性型)。球狀模型的幾何解釋，是因為它起源于2個半徑為a且球心距為2h的球體重疊部分的體積計算公式。它的一般的表示為0h=0；y(h)=曲+?非一；務(wù))，O〈hS〔Co+Cfh>a其中，二:為塊金常數(shù)，二:+C為基臺值，C為拱高，a為變程。當(dāng)。=0，C=1時，稱為標(biāo)準(zhǔn)球狀模型。其中，。為塊金常數(shù)，二:+C為基臺值，C為拱高，a為變程。(2)指數(shù)函數(shù)模型(在原點處為線性型)。它的一般公式為y(h)=/+c(i-廠")。此處a不是變程。因為當(dāng)h=3a時，有1-e-3^0.95^1。所以y(hR二:+C，故其變程為3a。當(dāng)二:=0，C=1時稱為標(biāo)準(zhǔn)指數(shù)函數(shù)模型。高斯模型(在原點處為拋物線型)。它的一般公式為V(h)=Co+C(1-eF)。此處a亦不是變程。因為當(dāng)h=-..3a時，有1-e-3^0.95^1，所以y(h)R二:+C，故其變程為3a。當(dāng)二:=0，C=1時稱為標(biāo)準(zhǔn)指數(shù)函數(shù)模型。6Kriging方法插值設(shè)Z(x)是點x承載的區(qū)域化變量，且是二階平穩(wěn)(或本征)的。Zi(i=1，2,…，〃)是一組離散的信息樣品數(shù)據(jù)，即重金屬濃度，它們是定義在點承載xi(i=1，2，…，n)上的?，F(xiàn)要對點x0承載處的區(qū)域化變量進行估計，所用的估計量為、二—.=，它是n個數(shù)值的加權(quán)線性組合。Kriging方法的原則，就是在保證這個估計量是無偏的，且估計方差最小的前提下，求出n個權(quán)值系數(shù)。(1)無偏性條件若要使二為二的無偏估計量，即要求E［二-二］=0，由此得到無偏條件：(2)Kriging方程組區(qū)域化變量在滿足二階平穩(wěn)的條件下推導(dǎo)，可以得到估計方差的計算公式

n+口i=l噸=g-z0]2-{g-z0])2=n+口i=l估計方差:三-對『?.的偏導(dǎo)數(shù)為寸=三三;二］;-2C:。在無偏性條件下，為了使估計方差最小，這是個求條件極值的問題，要用到拉格朗日乘子法。令。這里F是n個權(quán)系數(shù).".和的(n+1)元函數(shù)，求出F對.".(i=1，2，，n)和-?的偏導(dǎo)數(shù)，并令其為零，便得到下列Kriging方程組整理得整理得，，如果區(qū)域化變量只滿足本征假設(shè)，而不滿足二階平穩(wěn)假設(shè)，則利用協(xié)方差函數(shù)和變差函數(shù)的關(guān)系C(h)=C(0)-y(h)，可得用變差函數(shù)表示的Kriging方程組這里，yij=y(xIX)=y(x-Xj)。具體算法實現(xiàn)過程通過數(shù)學(xué)軟件matlab編程實現(xiàn)，(matlab代碼見附件)。⑶kriging模型的具體實現(xiàn)

隨機取100個觀測點作為原始數(shù)據(jù)，其分布如下圖所示:值得指出的是，在初步的數(shù)據(jù)處理中，可以發(fā)現(xiàn)該城區(qū)的地形并不是規(guī)則的長方形，但是為了方便數(shù)據(jù)處理，這里在一個30*20的區(qū)域內(nèi)隨機取點，剛好可以完全涵蓋該城區(qū)，但是這并不影響最終的結(jié)果，因為可以再得到模型之后再將不屬于城區(qū)的區(qū)域挖去，這在matlab軟件中是容易實現(xiàn)的。根據(jù)實驗變差函數(shù)結(jié)果，作數(shù)據(jù)擬合，選擇適當(dāng)?shù)淖儾詈瘮?shù)模型進行插值，確定變差函數(shù)計算公式。這里有3種常用變差函數(shù)模型：球狀模型、指數(shù)模型和高斯模型。但是這3種模型有3個關(guān)鍵參數(shù)需要確定，變程。、拱高C和塊金常數(shù)C0。本題中塊金常數(shù)C0=0，拱高C和變程a由實驗變差函數(shù)圖人工確定。對于第一種重金屬As，由實驗變差函數(shù)100個觀測點的觀測值得到的變差函數(shù)，可以計算得到"：.二；如圖：可以求得近似值：C=20和a=3。然后確定變差函數(shù)模型，選定變差函數(shù)模型后，變差函數(shù)的計算公式可以顯式寫出。這樣Kriging方程組的系數(shù)矩陣與增廣矩陣都已經(jīng)確定。然后進行Kriging插值，也就是求解Kriging方程組，確定加權(quán)系數(shù)，然后進行線性加權(quán)，即可得到該模型下的Kriging插值結(jié)果。（具體實現(xiàn)

見附件matlab代碼）經(jīng)畫圖比較，在本題中高斯變差函數(shù)模型的符合程度較好，建立高斯變差函數(shù)模型之后得到如下圖所示的函數(shù)圖像：1號金屬As分布原始數(shù)據(jù)結(jié)果1號金屬As分布原始數(shù)據(jù)結(jié)果1號金屬As分布kriging插值結(jié)果同樣的方式，對第二種重金屬Cd濃度做變差處理，值得指出的是，為方便比較，這里不再進行隨機取100個觀測點的工作，而是直接利用在計算第一種金屬As變差函數(shù)時候的隨機點數(shù)據(jù)，在此基礎(chǔ)上得到變差函數(shù)觀測數(shù)據(jù)"：.二；，如下圖：由實驗變差函數(shù)觀測圖可以近似求得C=126000和a=4。然后確定變差函數(shù)模型，選定變差函數(shù)模型后，變差函數(shù)的計算公式可以顯式寫出。這樣Kriging方程組的系數(shù)矩陣與增廣矩陣都已經(jīng)確定。然后進行Kriging插值，也就是求解Kriging方程組，確定加權(quán)系數(shù)，然后進行線性加權(quán)，即可得到該模型下的Kriging插值結(jié)果。經(jīng)畫圖比較，在本題中高斯變差函數(shù)模型的符合程度較好，建立高斯變差函數(shù)模型之后得到如下圖所示的圖像：2號金屬分布原始數(shù)據(jù)2號金屬kriging插值數(shù)據(jù)相應(yīng)的，利用2號金屬分布原始數(shù)據(jù)2號金屬kriging插值數(shù)據(jù)m號金屬空間分布4號金屬空間分布m號金屬空間分布4號金屬空間分布7確定各重金屬評價指標(biāo)的權(quán)重系數(shù)在綜合評定該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度時，還要對各重金屬評價指標(biāo)分別加權(quán)。權(quán)重是衡量因子集中某一因子對土壤污染程度影響相對大小的量，權(quán)重系數(shù)越大，則該因子對土壤的影響程度越大，在這里利用熵權(quán)法來確定各重金屬評價指標(biāo)的權(quán)重系數(shù)，具體原理及操作如下：方法優(yōu)點：客觀賦權(quán)法。背景：設(shè)有m個待評方案，n項評價指標(biāo)，形成原始指標(biāo)數(shù)據(jù)矩陣X=*）財（這里m=319，n=8），對于某項指標(biāo)Xj，指標(biāo)值Xij的值越大，則該指標(biāo)在綜合評價中所起的作用。在信息論中信息熵e=-k^Plnpiijij表示系統(tǒng)的有序程度，一個系統(tǒng)的有序程度越高，則信息熵越大，反之，一個系統(tǒng)的無序程度越高，則信息熵越小。所以，可以根據(jù)各項指標(biāo)的指標(biāo)值的差異程度，利用信息熵的這個工具計算出各指標(biāo)的權(quán)重。數(shù)據(jù)處理［氣，七…，氣］.?'X,X,...,Xmlm2mn上（i=1,2,...,m;j=1,2,...,n,），這里取n=8，m=319。由于參與評價的各項指標(biāo)有越大越優(yōu)型、越小越優(yōu)型，故需對矩陣中的特征值進行歸一化處理，方法如下：｛工『=#訂加以㈤）越大越優(yōu)型=X"min（虬J越小越優(yōu)型據(jù)此得到歸一化矩陣X’:X「X「…,X」X廣X21'，X22"X2n'X',X',..,Xmlm2mn(i=1,2,...,m;j=1,2,...,n,)這里以越大越優(yōu)型進行計算求解。計算第z?個因素下第j個評價值的比重PijP=-j'并x?

j=1計算第i個因素的熵值eie=-k芝Plnp,j=1

若取k=「―，則0<eJ1計算第i個因素的差異系數(shù)gi對于給定的e,越大，因素評價值的差異性越小，則因素在綜合評價中所起的作用越小。定義差異系數(shù)gj=1-e,，則當(dāng)因素g,越大時，因素越重要。定義權(quán)數(shù)七=」￥—，則w就是熵權(quán)法確定的權(quán)重。'jj=1記As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的權(quán)重矩陣為W=[W],w2,w3,w4,w5,w6,w7,w8]計算得到：W=[0.0233350.0475360.0651130.1831340.490280.0222410.0445820.123778]得到各取樣點加權(quán)綜合濃度指標(biāo)D=XW=[d1，d2，.,d319],；進一步得到i類區(qū)加權(quán)綜合濃度指標(biāo)d=-犬dij=1，其中，dj為所有屬于i類區(qū)的d，ni為其采樣點個數(shù)。得到評價結(jié)果得到五個區(qū)的加權(quán)綜合濃度指標(biāo)如下:區(qū)12345加權(quán)綜合濃度指標(biāo)值0.0021670.0059420.0009530.0040960.001787對各重金屬濃度背景值及標(biāo)準(zhǔn)差做相同歸一化處理，得到歸一化之后的背景平均值及標(biāo)準(zhǔn)差，然后對8種金屬歸一化之后的背景平均值和標(biāo)準(zhǔn)差進行加權(quán)求和，求得和各點加權(quán)綜合濃度指標(biāo)二.相對應(yīng)的綜合背景平均值和標(biāo)準(zhǔn)差如下：平均值對應(yīng)指標(biāo)標(biāo)準(zhǔn)差對應(yīng)指標(biāo)0.0007990.000195在此基礎(chǔ)上得到對污染程度等級進行分類的值如下：設(shè)輕度污染標(biāo)準(zhǔn)=平均值對應(yīng)指標(biāo)+3*標(biāo)準(zhǔn)差對應(yīng)指標(biāo);中度污染標(biāo)準(zhǔn)=平均值對應(yīng)指標(biāo)+8*標(biāo)準(zhǔn)差對應(yīng)指標(biāo);重度污染標(biāo)準(zhǔn)=平均值對應(yīng)指標(biāo)+15*標(biāo)準(zhǔn)差對應(yīng)指標(biāo);即:污染程度無輕度污染中度污染重度污染得分<0.0013850.001385~0.0023610.002361~0.003727>0.003727故可得到各區(qū)的污染等級如下表:區(qū)12345生活區(qū)工業(yè)區(qū)山區(qū)主干道公園綠地污染等級輕度污染重度污染無污染重度污染輕度污染由題設(shè)條件可知，按照功能劃分，城區(qū)一般可分為生活區(qū)、工業(yè)區(qū)、山區(qū)、主干道路區(qū)及公園綠地區(qū)等，分別記為1類區(qū)、2類區(qū)、……、5類區(qū)，不同的區(qū)域環(huán)境受人類活動影響的程度不同，通過以上的計算，得出的結(jié)論為：工業(yè)區(qū)與主干道嚴(yán)重污染，生活區(qū)與公園綠地輕度污染，山區(qū)幾乎不受污染，這與實際情況符合的較好。根據(jù)該金屬元素空間分布模型，可以確定工業(yè)區(qū)污染較為嚴(yán)重，故可以初步確定重金屬污染的主演原因在于工業(yè)區(qū)活動對環(huán)境的影響。

1.81.61.41.20.80.60.4■：;O'-0i104采樣點分布圖寺△O+生活區(qū)工業(yè)區(qū)山區(qū)主干道路區(qū)公園綠地去O松'*2Q,O*0。AA****X△□OA.-.勿半****2OOOOAA°AOAaoa°oxo半AA+O：：?：：OO°o°AXC?c-AA1.81.61.41.20.80.60.4■：;O'-0i104采樣點分布圖寺△O+生活區(qū)工業(yè)區(qū)山區(qū)主干道路區(qū)公園綠地去O松'*2Q,O*0。AA****X△□OA.-.勿半****2OOOOAA°AOAaoa°oxo半AA+O：：?：：OO°o°AXC?c-AA**i::::0.51.52.5x10得到污染源分布如下:各金屬污染源分布121110號號號號號號號號12-346678。口工寺■(>:△*+O2ud2uF(x,y)+k3,y)^—)—f3,y,u)=0ox2oy2其中，F(xiàn)(x,y)為污染源強度，K(x,y)為(x,y)處土壤傳播污染物的特性，f(x,y,u)為耗散強度。建立離散模型：設(shè)污染源為點源F(x.,y.)只有在源體必)處等于污染物的排放強度，在源處，恒穴假設(shè)各分區(qū)內(nèi)的土壤特性相同。則在非源處，82u+d2u)dx28y2ij則各區(qū)的kj=-!-Zkni其中，kj為j類區(qū)的平均土壤特性，%為i采樣點處的土壤特性。因為城區(qū)內(nèi)土壤特性相同，故f(x,y,u)可以簡寫為f(u)。取濃度場上的節(jié)點，建立離散方程組，通過求解方程組，可以分別得到F(x,y)，f(u)的一組數(shù)值，進行擬合，得出解析式，將解析式帶回擴散方程，即可得出擴散規(guī)律。六、模型評價與改進1模型評價為了給出8種主要重金屬的空間分布，并分析城區(qū)不同區(qū)域的污染程度。首先對空間位置比較散亂的采樣點此進行網(wǎng)格化的數(shù)據(jù)處理，然后通過Kriging方法進行空間散亂點的插值處理，在空間數(shù)據(jù)場滿足給定的統(tǒng)計分布特征的前提假設(shè)下進行插值。在求得相關(guān)參數(shù)后建立顯示的高斯模型，通過對函數(shù)的插值結(jié)果的觀察，與原始采樣數(shù)據(jù)的空間分布相比較，可以發(fā)現(xiàn)二者具有較好的吻合度，結(jié)果令人滿意。對城區(qū)內(nèi)不同區(qū)域重金屬污染程度的綜合評價，要對各重金屬評價指標(biāo)分別加權(quán)。利用熵權(quán)法來確定各重金屬評價指標(biāo)的權(quán)重系數(shù)，熵權(quán)法的實際意義在這里體現(xiàn)得尤為明顯，根據(jù)熵權(quán)法得到的相關(guān)系數(shù)均為正值，這一點也驗證了熵權(quán)法在尋找個金屬污染物權(quán)重時的正確性，然后由綜合權(quán)重進行線性加和，得到各個區(qū)域的綜合評定指標(biāo)，同時根據(jù)金屬含量背景值進行等級標(biāo)準(zhǔn)的劃分，從而確定不同區(qū)域的污染程度，結(jié)果與實際完全符合，說明熵權(quán)法的運用是正確的，從而便于找到重金屬污染的主要原因。2改進方向在建立重金屬污染的傳播特征模型，先假設(shè)了污染源的位置，然后考慮根據(jù)擴散定律建立模型，根據(jù)一維擴散方程建立模型，但是這樣的話在數(shù)據(jù)處理上必須首先根據(jù)采樣點濃度特征大致確定污染源的位置，然后建立方程，根據(jù)采樣數(shù)據(jù)提取信息求解，由于采樣點較多，本問題的處理可能會遺漏部分有用信息。如果與實際偏差較大，沒有得到有效結(jié)果，可以進一步尋求其他解決方法。對于問題四，可以基于問題一二三得到的模型進行綜合考慮，從數(shù)據(jù)的處理方式，以及地質(zhì)演變過程中時間變量的影響等因素入手，尋求更好的建立模型的方式。若引入時間變量，將擴散模型修改為：F3,J）+k3,j）（竺+孕）—f3,j,u）=4ox2dy2dt如此一來，還需要知道濃度的時間變化率，即不同時刻采樣點的濃度數(shù)值。然后要考慮二維（或三維）變差函數(shù)。但它們都是以一維函數(shù)為基礎(chǔ)的，只不過要考慮各向同性或各向異性，以及結(jié)構(gòu)的套合?；蛘呖梢赃M一步的加強考慮海拔Z的影響，得到四元的方程模型，然后建立求解。參考文獻周曉云，朱心雄.散亂數(shù)據(jù)點三角剖分方法綜述[J].工程圖學(xué)學(xué)報，1993,(01).王靖波，潘懋，張緒定.基于Kriging方法的空間散亂點插值[J].計算機輔助設(shè)計與圖形學(xué)學(xué)報，1999,(06).喬家君.改進的熵值法在河南省可持續(xù)發(fā)展能力評估中的應(yīng)用[J].資源科學(xué)，2004,(01).附件附件一重金屬空間分布模型的建立與求解的matlab實現(xiàn)注：為減少篇幅，此處代碼僅僅實現(xiàn)了對第一種金屬元素含量的空間分布的求解，其余個元素的求解方式類似可以得到。%坐標(biāo)網(wǎng)格化%zuobiao為采樣點坐標(biāo)矩陣%xx將坐標(biāo)單位化為公里x=zuobiao(:,1)';y=zuobiao(:,2)';z=zuobiao(:,3)';xx=x/1000;yy=y/1000;[X,Y]=meshgrid(0:0.5:29,0:0.5:19);Z=griddata(xx,yy,z,X,Y'v4');surf(X,YZ);shadinginterp%金屬濃度數(shù)據(jù)網(wǎng)格化jinshu_1=jinshu(:,1)';J_1=griddata(xx,yy,jinshu_1,X,Y'v4');jinshu_2=jinshu(:,2)';J_2=griddata(xx,yy,jinshu_2,X,Y'v4');jinshu_3=jinshu(:,3)';J_3=griddata(xx,yy,jinshu_3,X,Y'v4');jinshu_4=jinshu(:,4)';J_4=griddata(xx,yy,jinshu_4,X,Y'v4');jinshu_5=jinshu(:,5)';J_5=griddata(xx,yy,jinshu_5,X,Y'v4');jinshu_6=jinshu(:,6)';J_6=griddata(xx,yy,jinshu_6,X,Y'v4');jinshu_7=jinshu(:,7)';J_7=griddata(xx,yy,jinshu_7,X,Y'v4');jinshu_8=jinshu(:,8)';J_8=griddata(xx,yy,jinshu_8,X,Y'v4');%注：產(chǎn)生100個隨機坐標(biāo)的代碼只在對第一號金屬進行變差函數(shù)觀測的時候運行，在對其%他金屬進行變差函數(shù)觀測求解的時候只需利用當(dāng)前的100個觀測點即可，該代碼不再運行%產(chǎn)生100個隨機坐標(biāo)%s_suiji即為隨機觀測點得坐標(biāo)%fori=1:100%x_suiji(i,1:2)=[ceil(rand*30),ceil(rand*20)];%end%求變差函數(shù)的原始數(shù)據(jù)%bianchahanshu為變差函數(shù)矩陣bianchahanshu=[0,0];forj=1:100fork=1:100bianchahanshu=[bianchahanshu;...[sqrt((x_suiji(j,1)-x_suiji(k,1)).A2+(x_suiji(j,2)-x_suiji(k,2)).A2),...(J_1(x_suiji(j,2),x_suiji(j,1))-J_1(x_suiji(k,2),x_suiji(k,1))).A2]];endend%對變差函數(shù)矩陣進行處理，得到最終的變差函數(shù)矩陣form=1:600flag=find(bianchahanshu(:,1)==bianchahanshu(m,1));bianchahanshu(m,:)=sum(bianchahanshu(flag(:,1),:))/sum(flag==flag);bianchahanshu(flag(2:end,1),:)=[];mendbianchahanshu(:,2)=bianchahanshu(:,2)/2;%做出變差函數(shù)觀測點的圖像scatter(bianchahanshu(:,1),bianchahanshu(:,2))%根據(jù)變差函數(shù)，求每一個估計點權(quán)系數(shù)的系數(shù)矩陣以及增廣矩陣%C為系數(shù)矩陣%C_plus為增廣矩陣C=zeros(100,100);forp=1:100forq=1:100h=sqrt((x_suiji(p,1)-x_suiji(q,1)).A2+(x_suiji(p,2)-x_suiji(q,2)).A2);C(p,q)=20*(1-exp(-h/3));endendC_p