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

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

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

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

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

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

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

1.81.61.41.20.80.60.4■：;O'-0i104采樣點分布圖寺△O+生活區工業區山區主干道路區公園綠地去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+生活區工業區山區主干道路區公園綠地去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(x,y)為污染源強度，K(x,y)為(x,y)處土壤傳播污染物的特性，f(x,y,u)為耗散強度。建立離散模型：設污染源為點源F(x.,y.)只有在源體必)處等于污染物的排放強度，在源處，恒穴假設各分區內的土壤特性相同。則在非源處，82u+d2u)dx28y2ij則各區的kj=-!-Zkni其中，kj為j類區的平均土壤特性，%為i采樣點處的土壤特性。因為城區內土壤特性相同，故f(x,y,u)可以簡寫為f(u)。取濃度場上的節點，建立離散方程組，通過求解方程組，可以分別得到F(x,y)，f(u)的一組數值，進行擬合，得出解析式，將解析式帶回擴散方程，即可得出擴散規律。六、模型評價與改進1模型評價為了給出8種主要重金屬的空間分布，并分析城區不同區域的污染程度。首先對空間位置比較散亂的采樣點此進行網格化的數據處理，然后通過Kriging方法進行空間散亂點的插值處理，在空間數據場滿足給定的統計分布特征的前提假設下進行插值。在求得相關參數后建立顯示的高斯模型，通過對函數的插值結果的觀察，與原始采樣數據的空間分布相比較，可以發現二者具有較好的吻合度，結果令人滿意。對城區內不同區域重金屬污染程度的綜合評價，要對各重金屬評價指標分別加權。利用熵權法來確定各重金屬評價指標的權重系數，熵權法的實際意義在這里體現得尤為明顯，根據熵權法得到的相關系數均為正值，這一點也驗證了熵權法在尋找個金屬污染物權重時的正確性，然后由綜合權重進行線性加和，得到各個區域的綜合評定指標，同時根據金屬含量背景值進行等級標準的劃分，從而確定不同區域的污染程度，結果與實際完全符合，說明熵權法的運用是正確的，從而便于找到重金屬污染的主要原因。2改進方向在建立重金屬污染的傳播特征模型，先假設了污染源的位置，然后考慮根據擴散定律建立模型，根據一維擴散方程建立模型，但是這樣的話在數據處理上必須首先根據采樣點濃度特征大致確定污染源的位置，然后建立方程，根據采樣數據提取信息求解，由于采樣點較多，本問題的處理可能會遺漏部分有用信息。如果與實際偏差較大，沒有得到有效結果，可以進一步尋求其他解決方法。對于問題四，可以基于問題一二三得到的模型進行綜合考慮，從數據的處理方式，以及地質演變過程中時間變量的影響等因素入手，尋求更好的建立模型的方式。若引入時間變量，將擴散模型修改為：F3,J）+k3,j）（竺+孕）—f3,j,u）=4ox2dy2dt如此一來，還需要知道濃度的時間變化率，即不同時刻采樣點的濃度數值。然后要考慮二維（或三維）變差函數。但它們都是以一維函數為基礎的，只不過要考慮各向同性或各向異性，以及結構的套合。或者可以進一步的加強考慮海拔Z的影響，得到四元的方程模型，然后建立求解。參考文獻周曉云，朱心雄.散亂數據點三角剖分方法綜述[J].工程圖學學報，1993,(01).王靖波，潘懋，張緒定.基于Kriging方法的空間散亂點插值[J].計算機輔助設計與圖形學學報，1999,(06).喬家君.改進的熵值法在河南省可持續發展能力評估中的應用[J].資源科學，2004,(01).附件附件一重金屬空間分布模型的建立與求解的matlab實現注：為減少篇幅，此處代碼僅僅實現了對第一種金屬元素含量的空間分布的求解，其余個元素的求解方式類似可以得到。%坐標網格化%zuobiao為采樣點坐標矩陣%xx將坐標單位化為公里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%金屬濃度數據網格化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');%注：產生100個隨機坐標的代碼只在對第一號金屬進行變差函數觀測的時候運行，在對其%他金屬進行變差函數觀測求解的時候只需利用當前的100個觀測點即可，該代碼不再運行%產生100個隨機坐標%s_suiji即為隨機觀測點得坐標%fori=1:100%x_suiji(i,1:2)=[ceil(rand*30),ceil(rand*20)];%end%求變差函數的原始數據%bianchahanshu為變差函數矩陣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%對變差函數矩陣進行處理，得到最終的變差函數矩陣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;%做出變差函數觀測點的圖像scatter(bianchahanshu(:,1),bianchahanshu(:,2))%根據變差函數，求每一個估計點權系數的系數矩陣以及增廣矩陣%C為系數矩陣%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