高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽A題城市表層土壤重金屬污染分析

2011高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽承諾書(shū)我們仔細(xì)閱讀了中國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽的競(jìng)賽規(guī)則.我們完全明白，在競(jìng)賽開(kāi)始后參賽隊(duì)員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊(duì)外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問(wèn)題。我們知道，抄襲別人的成果是違反競(jìng)賽規(guī)則的,如果引用別人的成果或其他公開(kāi)的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻(xiàn)的表述方式在正文引用處和參考文獻(xiàn)中明確列出。我們鄭重承諾，嚴(yán)格遵守競(jìng)賽規(guī)則，以保證競(jìng)賽的公正、公平性。如有違反競(jìng)賽規(guī)則的行為，我們將受到嚴(yán)肅處理。我們參賽選擇的題號(hào)是（從A/B/C/D中選擇一項(xiàng)填寫(xiě)）：A 我們的參賽報(bào)名號(hào)為（如果賽區(qū)設(shè)置報(bào)名號(hào)的話）：所屬學(xué)校（請(qǐng)?zhí)顚?xiě)完整的全名）：浙江理工大學(xué)參賽隊(duì)員(打印并簽名)：1.馮世扣2.李漪3.陸天棟指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人(打印并簽名)：數(shù)模組日期：2011年9賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：2011高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽 編號(hào)專用頁(yè) 賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：賽區(qū)評(píng)閱記錄（可供賽區(qū)評(píng)閱時(shí)使用）：評(píng)閱人評(píng)分備注全國(guó)統(tǒng)一編號(hào)（由賽區(qū)組委會(huì)送交全國(guó)前編號(hào)）：全國(guó)評(píng)閱編號(hào)（由全國(guó)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：城市表層土壤重金屬污染分析摘要隨著城市經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷增加，人類活動(dòng)對(duì)城市環(huán)境質(zhì)量的影響日顯突出，因此研究城市地質(zhì)環(huán)境的變化十分重要。對(duì)于問(wèn)題一，我們首先對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，然后利用地統(tǒng)計(jì)學(xué)軟件surfer運(yùn)用反距離加權(quán)插值法（InverseDistancetoaPower）做出8種重金屬的空間等值線分布圖，緊接著我們用C語(yǔ)言編程算出五個(gè)功能區(qū)中各種重金屬含量的平均值。在分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度時(shí)，我們利用已經(jīng)統(tǒng)計(jì)的結(jié)果采用了單因子指數(shù)法，內(nèi)梅羅綜合污染指數(shù)法以及均方根綜合污染指數(shù)法進(jìn)行評(píng)價(jià)。用單因子評(píng)價(jià)方法得出了，該城區(qū)和Hg污染較為嚴(yán)重。用內(nèi)梅羅綜合污染指數(shù)法以及均方根綜合污染指數(shù)法得出五個(gè)功能區(qū)污染程度依次為工業(yè)區(qū)、主干道路區(qū)、生活區(qū)、公園綠地區(qū)和山區(qū)，其中工業(yè)區(qū)和主干道路區(qū)重度污染，生活區(qū)屬于中度污染，公園綠地區(qū)中度污染，山區(qū)為輕微污染。對(duì)問(wèn)題二，根據(jù)計(jì)算出來(lái)的單因子數(shù)值以及8張重金屬空間等值線分布圖，發(fā)現(xiàn)圖中有較多的地方含有高濃度的，且這些區(qū)域較為分散。這些區(qū)域可能是主干道路區(qū)，主要是因?yàn)槠囄矚庵泻写罅康摹Ｆ渌?張重金屬空間等值線分布圖分別只有幾個(gè)的區(qū)域含有高濃度的，Hg，，，，，，且這些區(qū)域集中是在工業(yè)區(qū)。因?yàn)楣I(yè)區(qū)排出來(lái)的“三廢”經(jīng)過(guò)大氣和水體最終污染土壤，使土壤中重金屬的濃度增加。問(wèn)題三我們把主成份分析法和擴(kuò)散問(wèn)題的偏微分方程結(jié)合來(lái)來(lái)建立模型。我們先用主成份分析法得到各重金屬間的相關(guān)性，并畫(huà)出主因子的空間等值線分布圖，從而先確定污染源的大概位置。然后再建立擴(kuò)散問(wèn)題的拋物型偏微分方程模型，為了確定源位置，通過(guò)微分方程的解析解和測(cè)量數(shù)據(jù)的關(guān)系，構(gòu)造了一個(gè)目標(biāo)函數(shù)優(yōu)化問(wèn)題，通過(guò)使用源附近濃度最大的點(diǎn)作為基準(zhǔn)點(diǎn)，巧妙地把非線性問(wèn)題轉(zhuǎn)化為線性問(wèn)題，用線性回歸確定了解析解的參數(shù)，并用SPSS軟件求解污染源的準(zhǔn)確位置。得出和來(lái)自同一個(gè)污染源（3863，7725，36），和來(lái)自共同污染源（2029，2717，34）（21524，11778，47），的污染源在（2016，5221，18）。的三個(gè)污染源分別在（6559，6947，34）（1356，5201，54）（17770，7515，54），的污染源（11209，9257，17）Hg污染源在（2894，2345，21）。對(duì)于問(wèn)題四由于擴(kuò)散問(wèn)題的偏微分方程模型是一個(gè)靜態(tài)的模型，因?yàn)榈叵滤淖饔?，重金屬在土壤里的傳播是一個(gè)流動(dòng)的模型。這時(shí)就要引入對(duì)流項(xiàng)，此時(shí)就會(huì)變成一個(gè)非常系數(shù)的偏微分方程，這種方程求解十分的復(fù)雜，可以采用反演的算法來(lái)求解。模型的創(chuàng)新性在于引入一個(gè)基準(zhǔn)點(diǎn)，把非線性的偏微分方程轉(zhuǎn)化為線性的偏微分方程。但是此模型是通過(guò)SPSS內(nèi)部的擬合的函數(shù)求得污染源位置，因此會(huì)存在一定的誤差。關(guān)鍵詞：重金屬土壤污染內(nèi)梅羅法主成分分析法擴(kuò)散問(wèn)題的偏微分方程一問(wèn)題重述1.1基本情況為了研究城市表層土壤重金屬污染程度，現(xiàn)對(duì)某城區(qū)地質(zhì)環(huán)境進(jìn)行勘察調(diào)查。按照功能劃分，城區(qū)一般可分為生活區(qū)、工業(yè)區(qū)、山區(qū)、主干道路區(qū)及公園綠地區(qū)等，分別記為1類區(qū)、2類區(qū)、……、5類區(qū)。我們將城區(qū)劃分為間距1公里左右的網(wǎng)格子區(qū)域，按照每平方公里1個(gè)采樣點(diǎn)對(duì)表層土（0～10厘米深度）進(jìn)行取樣、編號(hào)，并用GPS記錄采樣點(diǎn)的位置。應(yīng)用專門儀器測(cè)試分析，獲得了每個(gè)樣本所含的多種化學(xué)元素的濃度數(shù)據(jù)。另一方面，按照2公里的間距在那些遠(yuǎn)離人群及工業(yè)活動(dòng)的自然區(qū)取樣，將其作為該城區(qū)表層土壤中元素的背景值。1.2需要解決的問(wèn)題(1)給出8種主要重金屬元素在該城區(qū)的空間分布，并分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。(2)通過(guò)數(shù)據(jù)分析，說(shuō)明重金屬污染的主要原因(3)分析重金屬污染物的傳播特征，由此建立模型，確定污染源的位置。(4)分析你所建立模型的優(yōu)缺點(diǎn)，為更好地研究城市地質(zhì)環(huán)境的演變模式，還應(yīng)收集什么信息？有了這些信息，如何建立模型解決問(wèn)題？二符號(hào)說(shuō)明及基本假設(shè)2.1符號(hào)說(shuō)明（1）是指土壤中污染物的環(huán)境質(zhì)量指數(shù)；（2）是指污染物的實(shí)測(cè)值，單位為；（3）是指污染物的平均值，單位為；（4）則為采樣點(diǎn)處的濃度，為點(diǎn)源G的質(zhì)量；（5），，，為重金屬擴(kuò)散速度的算數(shù)平方根；（6）k為植物等對(duì)重金屬衰減系數(shù)的算術(shù)平方根；（7）表示采樣點(diǎn)的位置；（8）表示擬合出來(lái)的線性方程2.2基本假設(shè)（1）平穩(wěn)性假設(shè)，即假設(shè)在每個(gè)采樣點(diǎn)區(qū)域范圍內(nèi)重金屬污染物的濃度是均勻的，不考慮因某種特殊原因下某一點(diǎn)的某種重金屬含量顯著高于周圍的情況，如土壤中埋入了一節(jié)廢舊電池等。（2）假設(shè)同一種重金屬在不同的功能區(qū)的土壤傳播特征是一樣的（3）不考慮其他重金屬對(duì)該城區(qū)的污染（4）不考慮8種重金屬之間的化學(xué)物理反應(yīng)（5）假設(shè)附件中給出的數(shù)據(jù)都是準(zhǔn)確的三問(wèn)題分析3.1問(wèn)題1的分析為了給出8種主要重金屬元素在該城區(qū)的空間分布并且分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。題目中附錄一和附錄二提供了大量的數(shù)據(jù)，我們的任務(wù)是從這些龐大的調(diào)查數(shù)據(jù)中找出重金屬元素在不同區(qū)域的變化規(guī)律，這可以歸結(jié)為一個(gè)數(shù)理統(tǒng)計(jì)的問(wèn)題。只要我們把各個(gè)重金屬的位置和測(cè)量值統(tǒng)一標(biāo)在同一幅圖上，就可以得出重金屬在該城區(qū)的空間分布。對(duì)于分析不同城區(qū)重金屬的污染程度，我們引進(jìn)了新的評(píng)判指標(biāo)——均方根綜合污染指數(shù)法和內(nèi)梅羅綜合污染指數(shù)法，并且用編程求出城區(qū)內(nèi)不同區(qū)域的均方根綜合污染指數(shù)法和內(nèi)梅羅綜合污染指數(shù)，然后根據(jù)各自方法的評(píng)判指標(biāo)，得出該城區(qū)內(nèi)不同功能區(qū)重金屬的污染程度。3.2問(wèn)題2的分析對(duì)于問(wèn)題2，要分析重金屬的污染的主要原因，我們要結(jié)合重金屬在該城區(qū)的空間分布以及內(nèi)梅羅綜合污染指數(shù)法和均方根綜合污染指數(shù)法的結(jié)果，再查閱這幾年重金屬污染的相關(guān)文獻(xiàn)，找到了污染的主要原因3.3問(wèn)題3的分析要根據(jù)附錄中數(shù)據(jù)，確定污染源的位置。重金屬在土壤中的傳播屬于一個(gè)典型的擴(kuò)散問(wèn)題，但是由于該城區(qū)重金屬的采樣點(diǎn)較為的多，給我們用擴(kuò)散問(wèn)題的偏微分方程模型來(lái)尋找污染源造成了很大的困難。因此我們需要縮小污染源的范圍。另外由于不同重金屬之間可能存在相關(guān)性，即不同重金屬可能來(lái)自同一個(gè)污染源，因此我們想到了用主成份分析法。因?yàn)橹鞒煞莘治龇ú坏梢詼p少變量的個(gè)數(shù)，大大的減少了計(jì)算量，還可以畫(huà)出主因子的空間分布圖，這樣我們就可以找出重金屬較為集中的地方，也就是能夠?qū)さ轿廴驹吹拇蟾盼恢?。在確定污染源大概位置后，我們就可以用擴(kuò)散問(wèn)題的偏微分方程模型的方法來(lái)求出具體污染源的具體位置。四模型的建立與求解4.1基于內(nèi)梅羅綜合污染指數(shù)法的模型建立要建立8種重金屬在區(qū)域的空間分布，我們采用用反距離加權(quán)插值法（InverseDistancetoaPower）畫(huà)出每種重金屬的空間等值分布圖。在分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度過(guò)程中，我們先建立了內(nèi)梅羅綜合污染指數(shù)法評(píng)判模型。其中，單因子分指數(shù)的計(jì)算公式為：式中是指土壤中污染物的環(huán)境質(zhì)量指數(shù)；是指污染物的實(shí)測(cè)值，單位為；是指污染物的平均值，單位為。內(nèi)梅羅綜合污染指數(shù)的計(jì)算公式為：式中，為土壤污染中污染指數(shù)的最大值；是指土壤污染中污染指數(shù)的平均值。我們?cè)俳⒕礁C合污染指數(shù)法，其公式為：式中，P表示污染綜合指數(shù)；表示參加評(píng)價(jià)的污染物種類；表示污染物的實(shí)測(cè)值，；污染物的評(píng)價(jià)標(biāo)準(zhǔn)，；利用上述的公式算出各自的內(nèi)梅羅綜合指數(shù)和均方根指數(shù)后，再根據(jù)污染劃分指標(biāo)，分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。4.2內(nèi)梅羅綜合污染指數(shù)法評(píng)判模型求解4.2.1對(duì)8種主要重金屬元素在該城區(qū)的空間分布情況的分析圖4-1的空間等值線分布圖圖4-2的空間等值線分布圖從圖中可以看出，只有五個(gè)區(qū)域的濃度較為集中，其他地方的濃度都是較低的。因此可能是這個(gè)區(qū)域是加工工業(yè)區(qū)。對(duì)于表，圖中的有較多的地方含有高濃度的，且這些區(qū)域較為分散，這些可能是主干道路區(qū)。圖4-3的空間等值線分布圖圖4-4的空間等值線分布圖從圖中可以看出，只有兩個(gè)區(qū)域和的濃度較為集中，其他地方的濃度都是較低的。因此可能是這個(gè)區(qū)域是加工和工業(yè)區(qū)圖4-5Hg的空間等值線分布圖圖4-6的空間等值線分布圖從圖中可以看出，只有三個(gè)區(qū)域Hg的濃度較為集中，其他地方Hg的濃度都是較低的。因此可能是這個(gè)區(qū)域是Hg加工工業(yè)區(qū)。對(duì)于表，圖中的只有兩個(gè)區(qū)域的濃度較為集中，其他地方的濃度都是較低的，因此可能是這個(gè)區(qū)域是加工工業(yè)區(qū)圖4-7的空間等值線分布圖圖4-8的空間等值線分布圖對(duì)于表，圖中和的有較多的地方含有高濃度的，且這些區(qū)域較為分散，這些區(qū)域可能是主干道路區(qū)。4.2.2內(nèi)梅羅法和均方根法分析污染程度在模型的求解過(guò)程中，我們先將五個(gè)功能區(qū)的各種重金屬的數(shù)值帶入單因子分指數(shù)的計(jì)算公式：和內(nèi)梅羅公式：算出各個(gè)重金屬在不同功能區(qū)的單因子指數(shù)值和五個(gè)功能區(qū)的內(nèi)梅羅指數(shù)，數(shù)據(jù)見(jiàn)下表：表4-1各種金屬污染物在不同功能區(qū)的單因子指數(shù)值A(chǔ)sCdCrCuHgNiPbZn內(nèi)梅羅指數(shù)生活區(qū)1.742.661.492.233.433.17工業(yè)區(qū)2.003.021.729.6618.351.613.004.0313.533山區(qū)61.311.171.261.181.061.253主干道路區(qū)1.592.771.874.7112.771.432.053.529.427公園綠地區(qū)1.742.161.412.293.291.241.962.242.734表4-2土壤單因子指數(shù)分級(jí)標(biāo)準(zhǔn)污染等級(jí)未污染輕度污染中度污染重度污染極度污染研究單因子指數(shù)值可以發(fā)現(xiàn)，山區(qū)的污染最小，所有重金屬的單因子指數(shù)均在1和2之間，都屬于輕度污染。研究地土壤污染以Hg污染最為嚴(yán)重，其中在工業(yè)區(qū)和主干道路區(qū)Hg的單因子指數(shù)分別達(dá)到了18.35和12.77，都是屬于重度污染。生活區(qū)和公園綠地區(qū)分別則屬于中度污染，只有山地是輕度污染。五個(gè)功能區(qū)的和的單因子數(shù)值均在1到2之間，都是屬于輕度污染。五個(gè)功能區(qū)的除了生活區(qū)是中度污染外，其他地區(qū)均小于2，都屬于輕度污染。工業(yè)區(qū)的的單因子指數(shù)為9.66，已經(jīng)屬于極度污染，主干道路區(qū)屬于重度污染，生活區(qū)和公園綠地區(qū)為中度污染。工業(yè)區(qū)的的單因子指數(shù)為4.03,屬于重度污染，生活區(qū)，主干道路區(qū)和公寓綠地區(qū)均為中度污染。表4-3土壤綜合污染指數(shù)分級(jí)標(biāo)準(zhǔn)污染等級(jí)污染指數(shù)污染等級(jí)污染水平1安全清潔2警戒線尚清潔3輕污染土壤輕污染農(nóng)作物開(kāi)始受到污染4中污染土壤作物均受中度污染5重污染土壤作物均受污染且相當(dāng)嚴(yán)重通過(guò)內(nèi)梅羅指數(shù)我們可以很清楚的看到，工業(yè)區(qū)的污染是最嚴(yán)重的，內(nèi)梅羅指數(shù)為13.533，屬于重污染。主干道路區(qū)內(nèi)梅羅指數(shù)為9.427，也是屬于重污染。生活區(qū)的內(nèi)梅羅指數(shù)為3.170，也屬于重污染。公園綠地區(qū)的內(nèi)梅羅指數(shù)為2.734，屬于中污染。山地的污染最小，內(nèi)梅羅指數(shù)為1.253，屬于輕污染。然后我們將五個(gè)功能區(qū)的各種重金屬的數(shù)值帶入均方根公式：得出五個(gè)功能區(qū)的均方根指數(shù)，如下表:表4-4各功能區(qū)的均方根指數(shù)表功能區(qū)生活區(qū)工業(yè)區(qū)山區(qū)主干道路區(qū)公園綠地區(qū)均方根指數(shù)2.5747.6981.1945.2142.123表4-5土壤綜合污染指數(shù)分級(jí)標(biāo)準(zhǔn)污染等級(jí)污染指數(shù)污染等級(jí)污染水平1未污染清潔2輕微污染尚清潔3偏中度污染土壤輕污染農(nóng)作物開(kāi)始受到污染4中度污染土壤作物均受中度污染5偏重度污染土壤作物均受污染且相當(dāng)嚴(yán)重6重污染土壤作物均受污染且相當(dāng)嚴(yán)重7極度污染土壤作物均受污染且相當(dāng)嚴(yán)重通過(guò)均方根指數(shù)我們可以很清楚的看到，工業(yè)區(qū)的污染是最嚴(yán)重的，均方根指數(shù)為9.698，屬于極度污染。主干道路區(qū)均方根指數(shù)為5.214，也是屬于重污染。生活區(qū)的均方根指數(shù)為2.574，也屬于偏重度污染。公園綠地區(qū)的均方根指數(shù)為2.734，屬于偏重度污染。山地的污染最小，均方根指數(shù)為1.194，屬于輕微污染。4.3分析重金屬污染的主要原因土壤中的重金屬主要來(lái)源于成土母質(zhì)和人類活動(dòng)[1]。根據(jù)題目給出的土壤各重金屬的背景值和各個(gè)監(jiān)測(cè)點(diǎn)的實(shí)測(cè)值，再結(jié)合我們計(jì)算出來(lái)的單因子數(shù)值，我們可以判定，人類的活動(dòng)是造成這些功能區(qū)重金屬超標(biāo)的主要原因。這五個(gè)功能區(qū)污染程度依次為:工業(yè)區(qū)﹥主干道路區(qū)﹥生活區(qū)﹥公園綠地區(qū)﹥山區(qū)。工業(yè)區(qū)中，Hg，，，已經(jīng)嚴(yán)重超標(biāo)，另外三種也超標(biāo)，主要由于工業(yè)區(qū)有較多的金屬冶煉廠，蓄電池加工廠，印刷廠和化工廠等，這些工廠都需要大量的金屬原材料，煤炭和石油等，這些都會(huì)產(chǎn)生大量的重金屬?gòu)U渣廢氣。這些廢渣廢氣排到土壤，河流和大氣里，會(huì)造成工業(yè)區(qū)土壤重金屬含量嚴(yán)重超標(biāo)，另外這些廢氣廢水會(huì)隨著大氣水體影響附近的地區(qū)的土壤。比如蓄電池加工廠周圍的土壤中的含量會(huì)明顯超標(biāo)。不同的金屬冶煉廠對(duì)各重金屬污染在周邊土壤中的積累貢獻(xiàn)各異，如冶煉廠附近土壤中不僅的含量會(huì)增加，同時(shí)的含量也會(huì)有所增加[4]。主干道路區(qū)的,,,污染也較為嚴(yán)重，主要是由于機(jī)動(dòng)車會(huì)產(chǎn)生大量的尾氣。汽車尾氣中含有大量的，，，等重金屬元素，這些重金屬會(huì)嚴(yán)重影響主干道路區(qū)的土壤。通過(guò)對(duì)汽車尾氣顆粒物中重金屬含量的測(cè)定發(fā)現(xiàn)，的含量為37%，，，的含量分別為34.5%，22.6%，2.6%[2]。土壤中的元素是不能隨著時(shí)間被土壤降解的，這就導(dǎo)致元素會(huì)隨著時(shí)間慢慢積累。另外元素是汽車硬度的添加劑，在汽車輪胎摩擦產(chǎn)生的粉塵會(huì)引起土壤，含量的增多。生活區(qū)的污染主要是由于生活垃圾和生活污水的排放以及生活用煤。垃圾堆積在土地上，垃圾中的重金屬會(huì)滲透到土壤里，這就會(huì)導(dǎo)致土壤重金屬含量的增加。生活污水的排放，重金屬就會(huì)流入河流，最終影響生活區(qū)的土壤和其他地區(qū)的土壤。居民用煤也是產(chǎn)生重金屬的重要原因。煤炭燃燒過(guò)程中產(chǎn)生的易于富集在干灰中，且和更易于富集于細(xì)灰粒中，然后隨著細(xì)灰粒排放到大氣中，這些重金屬就會(huì)進(jìn)入大氣，最終沉降到土壤中。公園綠地區(qū)的污染較輕，其中，Hg，都是中度超標(biāo)。這有可能是受工業(yè)區(qū)的廢氣和汽車尾氣影響。因?yàn)榇髿庵械闹亟饘俳?jīng)過(guò)自然沉淀和雨淋會(huì)進(jìn)入土壤中，這就有可能影響到公園綠地區(qū)的土壤。另外旅客產(chǎn)生的垃圾也會(huì)影響公園綠地去的土壤。山區(qū)的污染在這幾個(gè)功能區(qū)里面是最輕的。山區(qū)的重金屬來(lái)源可能來(lái)自大氣，附近的河流和農(nóng)藥的使用。但是植物能夠降解重金屬，這會(huì)使得山區(qū)的土壤重金屬含量大大的降低。4.4基于擴(kuò)散問(wèn)題的偏微分方程建立首先我們先建立主成分分析法模型，來(lái)找出污染源的大概位置。第一步，我們先對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，把數(shù)據(jù)根據(jù)功能區(qū)重新排列，依次為1區(qū)，2區(qū)，3區(qū)，4區(qū)，5區(qū)。第二步我們?cè)倮密浖?jì)算出標(biāo)準(zhǔn)化數(shù)據(jù)的相關(guān)系數(shù)矩陣，并且算出相關(guān)系數(shù)矩陣的特征值，特征向量，方差貢獻(xiàn)率，累計(jì)方差貢獻(xiàn)率。正交的方差最大旋轉(zhuǎn)法來(lái)求解旋轉(zhuǎn)后的特征向量，方差貢獻(xiàn)率，累計(jì)方差貢獻(xiàn)率。通過(guò)旋轉(zhuǎn)后的累計(jì)方差貢獻(xiàn)率，我們就可以確定哪幾個(gè)因子為主因子，哪幾個(gè)因子為次因子。第三步我們?cè)龠M(jìn)行正交互換，使用方差最大法，求出旋轉(zhuǎn)因子的載荷矩陣。第四步我們通過(guò)分析旋轉(zhuǎn)因子的載荷矩陣，就可以知道哪些重金屬為主因子，哪些重金屬屬于次因子。為了更加精確的確定污染源的位置，我們做出了主因子在5個(gè)功能區(qū)的空間分布等直線圖，這樣能夠更加直觀的看出污染源大概在空間上的分布情況。然后我們就可以建立了擴(kuò)散問(wèn)題的偏微分方程模型，即（1）其中點(diǎn)G（，，）為假設(shè)的污染源擴(kuò)散點(diǎn)，則為G點(diǎn)處的濃度，，，，為重金屬擴(kuò)散速度的算數(shù)平方根，k為植物等對(duì)重金屬衰減系數(shù)的算術(shù)平方根，為點(diǎn)源G的質(zhì)量。（，，）為取樣點(diǎn)的位置。4.5擴(kuò)散問(wèn)題的偏微分模型的求解4.5在用主成份分析法的過(guò)程中，我們先用SPSS軟件得出了各種重金屬之間的相關(guān)性。相關(guān)性分析就是研究?jī)蓚€(gè)或兩個(gè)以上變量之間的相互關(guān)系的統(tǒng)計(jì)分析方法。在環(huán)境研究中常用來(lái)定性的分析研究各環(huán)境變量之間的相互關(guān)系，而變量之間的相互關(guān)系用相關(guān)系數(shù)和相關(guān)系數(shù)的顯著水平來(lái)表示[3]。如果相關(guān)系數(shù)為正，表示兩個(gè)變量正相關(guān)，取值越小相關(guān)性越強(qiáng)。相反，如果相關(guān)系數(shù)為負(fù)值，表示兩個(gè)變量負(fù)相關(guān)，值越小則表示相關(guān)性越強(qiáng)。土壤中重金屬含量的相關(guān)性可以推測(cè)出重金屬是不是來(lái)自同一個(gè)污染源[4]。如果重金屬含量有著很強(qiáng)的相關(guān)性，表明重金屬的來(lái)源可能相同，如果重金屬的相關(guān)性不是很強(qiáng)，則重金屬的來(lái)源可能不止一個(gè)。下表列出了各種重金屬之間的相關(guān)性：表4-6重金屬的相關(guān)性表AsCdCrCuHgNiPbZnAs10.2550.1890.160.0640.3170.290.247Cd0.25510.3520.3970.2650.3290.660.431Cr0.1890.35210.5320.1030.7160.3830.424Cu0.160.3970.53210.4170.4950.520.387Hg0.0640.2650.1030.41710.1030.2980.196Ni0.3170.3290.7160.4950.10310.3070.436Pb0.290.660.3830.520.2980.30710.494Zn0.2470.4310.4240.3870.1960.4360.4941從圖中我們可以看出各個(gè)重金屬之間的相關(guān)性不是很強(qiáng)。和的相關(guān)度為0.532，為中度極顯著正相關(guān)，和的相關(guān)度為0.716，為高度極顯著正相關(guān)，和均為中度極顯著正相關(guān)，和的相關(guān)度為0.660，為中度極顯著正相關(guān)，和的相關(guān)度為0.494，也為中度極顯著正相關(guān)。4.5由于該城區(qū)有八種重金屬的污染，而且有些金屬之間還存在一定的相關(guān)性，這就會(huì)導(dǎo)致放映的信息重疊。因此我們要降低因子數(shù)。而因子分析則是通過(guò)降維方法進(jìn)行簡(jiǎn)化得到綜合指標(biāo)，綜合指標(biāo)之間既互不相關(guān)，又能反映原來(lái)的觀察指標(biāo)的信息。我們用SPSS軟件算出了特征值，貢獻(xiàn)率和累計(jì)的貢獻(xiàn)率。但是主成份分析法只強(qiáng)調(diào)離差而不是強(qiáng)調(diào)變量在樣品中的百分含量，因此我們用正交的方差最大旋轉(zhuǎn)法使每一個(gè)主因子只與最少個(gè)數(shù)的變量有相關(guān)的關(guān)系[6]。算出來(lái)的旋轉(zhuǎn)變化后的特征值，貢獻(xiàn)率和累計(jì)貢獻(xiàn)率，如表所示：表4-7特征值和累計(jì)貢獻(xiàn)率表成份旋轉(zhuǎn)變化前旋轉(zhuǎn)變化后特征值貢獻(xiàn)率%累計(jì)貢獻(xiàn)率%特征值貢獻(xiàn)率%累積貢獻(xiàn)率%13.5644.544.51.82722.83422.83421.1514.37758.8771.52519.06341.89730.96512.06370.9411.04513.06854.96540.7689.59680.5371.03212.89867.86350.5787.2287.7561.0212.74480.60760.4325.39993.1561.00412.54993.156從表格里可以看出，累計(jì)貢獻(xiàn)率為93.156%的前提下，我們得到了6個(gè)主因子，也就是意味著6個(gè)主因子提供了93.156%的信息，滿足主成份分析法的條件。這樣大大地簡(jiǎn)化問(wèn)題，且能偶提供93.156%的信息，有較高的準(zhǔn)確性。另外我們發(fā)現(xiàn)在旋轉(zhuǎn)前后累計(jì)貢獻(xiàn)率都是為93.156%，表明總量信息沒(méi)有丟失，我們的模型具有較高的準(zhǔn)確性。

表方差極大正交旋轉(zhuǎn)前的因子載荷矩陣。對(duì)旋轉(zhuǎn)后的數(shù)據(jù)進(jìn)行更加相信的分析，我們發(fā)現(xiàn)，1和2的貢獻(xiàn)率大約在20%左右，這可以說(shuō)明這兩個(gè)因子可能是這個(gè)地區(qū)的最大的污染源，對(duì)這個(gè)地區(qū)的環(huán)境影響最大。而3，4，5，6四個(gè)因子的貢獻(xiàn)率大概在13%左右，表明這四個(gè)因子對(duì)這個(gè)地區(qū)的污染沒(méi)有1，2兩個(gè)因子那么重要，但是也是也對(duì)土壤有著較大的影響。表4-8旋轉(zhuǎn)前因子載荷矩陣123456As0.426-0.20.6810.551-0.026-0.065Cd0.7110.2810.282-0.322-0.2540.325Cr0.735-0.444-0.303-0.046-0.110.098Cu0.7560.125-0.3650.137-0.155-0.408Hg0.4080.673-0.2970.4490.1540.236Ni0.723-0.515-0.190.137-0.0140.2Pb0.7640.3140.237-0.248-0.158-0.217Zn0.699-0.0370.123-0.2410.654-0.06表4-9旋轉(zhuǎn)后的因子載荷矩陣123456As0.1350.130.0160.9740.0430.084Cd0.2230.9180.1480.0860.0310.126Cr0.8590.187-0.0160.0040.2450.141Cu0.3950.1770.2730.0310.810.11Hg0.0170.1370.9670.0160.1730.069Ni0.890.1010.0510.1960.1030.17Pb0.0730.7180.0610.1620.4960.275Zn0.2610.2350.0820.0970.1220.916結(jié)合上面兩幅表格一起分析可以得出：旋轉(zhuǎn)前后因子載荷矩陣數(shù)值較為接近。然后我們繼續(xù)分析旋轉(zhuǎn)后的因子載荷矩陣，發(fā)現(xiàn)重金屬與因子的載荷值越大，則表示重金屬和因子的關(guān)系越密切。例如和與因子1的載荷分別為0.859和0.890。我們就認(rèn)為因子1為和的組合。以此類推，因子2為和的組合，因子3為Hg，因子4為，因子5為，因子6為。因此和，和極有可能來(lái)自同一污染源。在前面的重金屬相關(guān)性分析中，和，和都具有較好的。另外的重金屬則來(lái)自不同的污染源。4.5.3首先我們運(yùn)用傅立葉變換出求出偏微分方程的通解，求出了的通解為（2）如果認(rèn)為經(jīng)過(guò)相當(dāng)長(zhǎng)的時(shí)間后，擴(kuò)散已經(jīng)終止，物質(zhì)分布處于平穩(wěn)狀態(tài)，則方程1中，于是方程1就可以有線性橢圓型方程的邊值問(wèn)題其中是區(qū)域D的邊界。由于點(diǎn)源G的質(zhì)量與位置，（，，）是已知的，前面我們用主成份分析法已經(jīng)找出來(lái)的大概污染源的位置。為了確定方程（2）中的7個(gè)參數(shù)，我們給出下面一個(gè)目標(biāo)函數(shù)：取樣本點(diǎn)，使得最小,考慮到解的表達(dá)式中含有指數(shù)取對(duì)數(shù)，得：其中包括了項(xiàng)，是非線性的。所以我們引入一個(gè)濃度最大的點(diǎn)作為基準(zhǔn)點(diǎn)，這樣可以消掉一個(gè)參數(shù)M，并把非線性問(wèn)題轉(zhuǎn)換為線性問(wèn)題，即求下面函數(shù)的最小值：現(xiàn)在我們先在大概污染源的中心取一點(diǎn)H，然后在大概污染源附近取一系列測(cè)量點(diǎn)：···。令t=1,把一系列的測(cè)量點(diǎn)和H點(diǎn)帶入偏微分方程的通解中，再把測(cè)量點(diǎn)的物質(zhì)密度與通解公式都取對(duì)數(shù)，得出兩條等式，兩條等式相減，則其中系數(shù)，，。我們把點(diǎn)源，點(diǎn)H，一些列的測(cè)量點(diǎn)···進(jìn)行多元一次線性方程擬合，可以擬合出一條依次線性方程。然后根據(jù)系數(shù)就可以求出，，的值，即確定了污染源的位置。表4-10污染源的位置坐標(biāo)位置AsCd和PbCr和NiCuHgZnx6559.792029.583863.542016289411209.1y6947.712717.967725.265221.1823459257.45z34.7384.1736.19182117.02x13560.5321524.8y5201.2111778.28z5447.64x17770.5y7515.97z54.284.6問(wèn)題4求解4.6.1（1）我們先用主成份分析法確定污染源的大概位置后，在用擴(kuò)散問(wèn)題的偏微分方程模型來(lái)求解，這樣大大的減少了建模過(guò)程中的計(jì)算量，并且能過(guò)很快的算出污染源的準(zhǔn)確位置。（2）此模型還考慮了重金屬在土壤中的衰減，比較符合實(shí)際，具有較好的有科學(xué)性。4.6.2（1）在建立模型三的過(guò)程中，我們考慮的是粉塵模型，沒(méi)有考慮水流的作用，這樣建立的數(shù)學(xué)模型與實(shí)際可能有一定的誤差。（2）由于本模型在求解過(guò)程中，通過(guò)SPSS軟件來(lái)擬合大量數(shù)據(jù)，從而求出污染源的位置。用擬合算法會(huì)存在一定的誤差。4.6.3在知道了模型三的缺點(diǎn)后，我們需要采集相關(guān)的信息加以完善。比如采集多個(gè)不同時(shí)刻重金屬在相同樣本點(diǎn)的濃度。在模型中加上對(duì)流系數(shù)后，模型就不再是重金屬在的靜止的土壤中傳播了，而是在相對(duì)流動(dòng)的土壤中傳播，這樣更加的符合實(shí)際，也能更加準(zhǔn)確的求出污染源的位置。模型三除了需要考慮土壤的對(duì)流系數(shù)，還可以考慮重金屬在空間的變異性和自流系數(shù)。我們考慮到地下水流的影響，此時(shí)原函數(shù)就變?yōu)椋何覀兩踔吝€應(yīng)該考慮到重金屬在各個(gè)方向上擴(kuò)散系數(shù)的變化，即a、b、c是關(guān)于時(shí)間的函數(shù)。此時(shí)該方程是個(gè)非常難解的非常系數(shù)的偏微分方程，很難把解析解求出來(lái)，已知樣本點(diǎn)的一些測(cè)量數(shù)據(jù)要求出源位置，這是偏微分方程反問(wèn)題研究的范疇。已知源項(xiàng)和邊值條件，由初邊值條件求解即是人們熟知的正問(wèn)題。而源項(xiàng)反問(wèn)題的提法是，已知在某些點(diǎn)的測(cè)量數(shù)據(jù)或某個(gè)時(shí)刻T的測(cè)量數(shù)據(jù)反求源項(xiàng)，且該反問(wèn)題在一般意義下是不適定的。由于該偏微分方程反問(wèn)題含有很多參數(shù)，直接得到解析解是很困難的，必須使用數(shù)值微分的方法做反演計(jì)算（見(jiàn)[8,9]）。五模型的結(jié)果分析及評(píng)價(jià)5.1內(nèi)梅羅綜合污染指數(shù)法評(píng)判模型的結(jié)果分析對(duì)于模型一，我們用內(nèi)梅羅綜合污染指數(shù)法和均方根綜合指數(shù)法來(lái)建立模型，雖然用內(nèi)梅羅法求出來(lái)的數(shù)值稍大于均方根法求出來(lái)的數(shù)值，但是整體上的污染等級(jí)還是較為相近的。用兩種方法建立模型求解的結(jié)果中，污染程度都是工業(yè)區(qū)>主干道路區(qū)>生活區(qū)>公園綠地區(qū)>山區(qū)。其中工業(yè)區(qū)和主干道區(qū)都是屬于重度污染！因此兩個(gè)模型都具有相當(dāng)好的可行性。5.2擴(kuò)散問(wèn)題的偏微分方程模型的結(jié)果分析我們用模型3求解出來(lái)的污染源的位置和第一問(wèn)中的重金屬的空間等值線分布圖進(jìn)行了比較，發(fā)現(xiàn)污染源附近重金屬濃度相當(dāng)?shù)母?。這也說(shuō)明模型三具有一定的準(zhǔn)確性。5.3模型的優(yōu)點(diǎn)（1）內(nèi)梅羅綜合污染指數(shù)法操作簡(jiǎn)單，運(yùn)算也較為方便，只需要算出內(nèi)梅羅綜合污染指數(shù)，然后再跟評(píng)判標(biāo)準(zhǔn)進(jìn)行比較，就可以得出地區(qū)的污染程度。另外內(nèi)梅羅綜合污染指數(shù)在加權(quán)過(guò)程中避免了權(quán)系數(shù)中主觀因素的影響，兼顧了單因子污染指數(shù)的平均值和最高值。具有一定的可靠性。（2）模型三的優(yōu)點(diǎn)在問(wèn)題四中已有詳盡的說(shuō)明。5.4模型的缺點(diǎn)（1）內(nèi)梅羅綜合污染法在計(jì)算過(guò)程中，過(guò)分的看重了極大值對(duì)水體的污染，如果有一項(xiàng)指標(biāo)的值很大，而其他值較小，這樣也會(huì)使內(nèi)梅羅綜合指數(shù)偏大[7]。另外內(nèi)梅羅綜合污染法還忽略了重金屬的毒性的差別。因此得到的結(jié)果往往比較機(jī)械。（2）模型三的缺點(diǎn)在問(wèn)題四中已有詳盡的說(shuō)明。在改進(jìn)內(nèi)梅羅綜合污染法的過(guò)程中，我們還考慮了不同重金屬的毒性對(duì)的影響，得到改進(jìn)后的內(nèi)梅羅綜合污染指數(shù)法，其公式為：其中為權(quán)重值最大的污染因子的與的比值。表示第種污染因子的相關(guān)性比值。令，其中為第種污染因子的權(quán)重值，。改進(jìn)之后的內(nèi)梅羅綜合污染法一定程度上消除了極大值的影響。并提高了重金屬毒性指標(biāo)的影響，能夠更加客觀的評(píng)價(jià)地區(qū)的污染狀況。參考文獻(xiàn)[1]張嘉穎，大寶山礦區(qū)橫石河流域重金屬污染及其生態(tài)評(píng)價(jià)3/WFknowledgeServer_Mirror/D/Thesis_Y1476990.aspx，2011.9.9[2]張志紅，楊文敏，汽油車排除顆粒物化學(xué)組分分析。中國(guó)公共衛(wèi)生。2001，17（7）：623-624[3]林艷，基于地統(tǒng)計(jì)學(xué)與GIS的土壤重金屬污染評(píng)價(jià)與預(yù)測(cè)，/p-68758645.html，2011.9.10[4]周燕，崇明島土壤重金屬污染的空間分析及污染評(píng)價(jià)，/kcms/detail/detail.aspx?dbCode=CMFD&dbName=CMFD2009&FileName=2009187292.nh&filetitle=，2011.9.10[5]張乃明，李保國(guó)，胡克林，太原污灌區(qū)土壤重金屬與鹽分含量的空間分布特征【J】．環(huán)境科學(xué)學(xué)報(bào)，21(3)：349～353，2001[6]王雄軍，賴健青，魯艷紅，李德勝，周繼紅，王建武，基于因子分析法研究太原市土壤重金屬污染的主要來(lái)源，17（2）：671-676，2008[7]李亞松，張兆吉，費(fèi)宇紅，王昭，內(nèi)梅羅指數(shù)評(píng)價(jià)法的修正及其應(yīng)用，25（6）：48-50，2009[8]YAMAMOTOM.Conditionalstabilityindeterminationofforcetermsofheatequationsinarectangle[J].Math.Comput.Modelling,18:79-88,1993[9]CHOULLIM,YAMAMOTOM.Conditionalstabilityindeterminingaheatsource[J].JournalofInverseandIll-PosedProblems,12:233-336,2004附錄求解內(nèi)梅羅指數(shù)的程序#include<stdio.h>#include<math.h>intmain(){doublep[5][8];doublec[5][8]={6.27,289.96,69.02,49.4,93.04,18.34,69.11,237.01, 7.25,393.11,53.41,127.54,642.36,19.81,93.04,277.93, 4.04,152.32,38.96,17.32,40.96,15.45,36.56,73.29,5.71,360.01,58.05,62.21,446.82,17.62,63.53,242.85,6.26,280.54,43.64,30.19,114.99,15.29,60.71,154.24};doubles[8]={3.6,130,31,13.2,35,12.3,31,69};inti,j;for(i=0;i<5;i++) for(j=0;j<8;j++) p[i][j]=c[i][j]/s[j];doublesum,ave,max,temp;doubleans[5];for(i=0;i<5;i++){sum=p[i][0];max=p[i][0];for(j=1;j<8;j++) { if(p[i][j]>max) max=p[i][j]; sum=sum+p[i][j]; } ave=sum/8.0; ans[i]=sqrt((max*max+ave*ave)/2.0);}for(i=0;i<5;i++)printf("%lf\n",ans[i]);return0;}求解單因子指數(shù)的程序#include<stdio.h>#include<string.h>intmain(){doubles[319][9];doublesum[6][9];doubless[8]={3.6,130,31,13.2,35,12.3,31,69};intcount[6]={0,0,0,0,0,0};FILE*fp;fp=fopen("1.txt","r");inti,j;memset(sum,'0',sizeof(sum));for(i=0;i<319;i++){ fscanf(fp,"%lf%lf%lf%lf%lf%lf%lf%lf%lf",&s[i][0],&s[i][1],&