1、城市表層土壤重金屬污染分析摘要針對(duì)第一問(wèn)，本文利用 griddata 函數(shù)中的 cubic 插值，畫(huà)出整個(gè)城區(qū)內(nèi)各重金屬污染分布等值線圖，以最直觀的方式，大致確定重金屬元素的空間分布。通過(guò)計(jì)算采樣點(diǎn)對(duì)各金屬元素地質(zhì)累積指數(shù)，以及在此基礎(chǔ)上修正的量不同區(qū)域重金屬污染程度的評(píng)價(jià)體系。綜合指數(shù)評(píng)價(jià)法，建立起衡針對(duì)第二問(wèn)，對(duì)于一種金屬元素，通過(guò)比較不同功能區(qū)域的污染情況，并結(jié)合該金屬在現(xiàn)實(shí)中的實(shí)際應(yīng)用，確定該金屬元素的來(lái)源與原因（工業(yè)廢水、生活污水、汽車(chē)尾氣，固體廢棄物等）。第三問(wèn)為源項(xiàng)識(shí)別反問(wèn)題。本文先建立重金屬濃度擴(kuò)散的偏微分方程，并轉(zhuǎn)化成差分模型，再引入相對(duì)海拔對(duì)擴(kuò)散速度影響因子。從數(shù)據(jù)中充分挖

2、掘信息，估計(jì)出模型中的相關(guān)參量。在污染區(qū)域任選一點(diǎn)(x, y) 作為污染源，則問(wèn)題轉(zhuǎn)化成已知源和初值狀態(tài)的正問(wèn)題，求解并 到采樣點(diǎn)實(shí)際測(cè)量值與模型求解值之間的偏差 Dev 。以 minDev為目標(biāo)函數(shù)，讓(x, y) 遍歷所有可能解，從而確定出污染源的坐標(biāo)位置，并在 SPSS 中，用 Pearson 相關(guān)性分析檢驗(yàn)?zāi)Ｐ?。受限于題目數(shù)據(jù)匱乏，第三問(wèn)所建模型較為粗糙。第四問(wèn)中，若能通過(guò)取樣和實(shí)驗(yàn)測(cè)得相關(guān)參數(shù)，可以建立起一個(gè)涉及擴(kuò)散、對(duì)流和衰減的濃度變化偏微分方程，據(jù)此可求的更加精確的結(jié)果。差分方程 Pearson 相關(guān)性分析：重金屬反問(wèn)題11.問(wèn)題重述問(wèn)題背景：隨著城市經(jīng)濟(jì)的快速發(fā)展和城市的不斷增

3、加，人類(lèi)活動(dòng)對(duì)城市環(huán)境質(zhì)量的影響日益突出。其中，城市表層土壤重金屬污染也日趨顯著。對(duì)城市土壤地質(zhì)環(huán)境異常的查證，以及如何應(yīng)用查證獲得的海量數(shù)據(jù)資料開(kāi)展城市環(huán)境類(lèi)活動(dòng)影響下城市地質(zhì)環(huán)境的演變模式，日益成為人們關(guān)注的焦點(diǎn)。價(jià)，研究人城市分區(qū)：按照功能劃分，城區(qū)一般可分為生活區(qū)、工業(yè)區(qū)、山區(qū)、主干道路區(qū)（交通區(qū)）及公園綠地區(qū)，分別記為 1 類(lèi)區(qū)、2 類(lèi)區(qū)、類(lèi)活動(dòng)影響的程度不同。5 類(lèi)區(qū)，不同的區(qū)域環(huán)境受人土壤采樣：現(xiàn)對(duì)某城市城區(qū)土壤地質(zhì)環(huán)境。為此，將所的城區(qū)劃分為間距 1 公里左右的網(wǎng)格子區(qū)域，按照每平方公里 1 個(gè)采樣點(diǎn)對(duì)表層土（010 厘米深度），并用 GPS進(jìn)行取樣、采樣點(diǎn)的位置。應(yīng)用專(zhuān)門(mén)儀器

4、測(cè)試分析，獲得了每個(gè)樣本所含的多種化學(xué)元素的濃度數(shù)據(jù)。另一方面，對(duì)遠(yuǎn)離人群及工業(yè)區(qū)的自然區(qū)進(jìn)行取樣，并將其作為該城區(qū)表層土壤的背景值。為了了解城區(qū)主要重金屬的分布，對(duì)城區(qū)環(huán)境的質(zhì)量標(biāo)準(zhǔn)進(jìn)行評(píng)價(jià)，同時(shí)研究人類(lèi)活動(dòng)影響下城市地址環(huán)境的演變模式，需要完成以下四個(gè)任務(wù)：1)給出 8 中主要重金屬元素在該城區(qū)的空間分布，并分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度；通過(guò)數(shù)據(jù)分析，說(shuō)明重金屬污染的主要原因；2)3)4)分析重金屬污染物的特征，并建立污染物模型，確定污染源位置；評(píng)價(jià)所建立模型的優(yōu)缺點(diǎn)，考慮還應(yīng)收集哪些信息以及利用這些信息建立更完善的模型以深入研究城市地質(zhì)環(huán)境的演變。2.基本假設(shè)1)2)3)一個(gè)污染

5、區(qū)域只由一個(gè)污染源產(chǎn)生。重金屬污染只在表層土壤中擴(kuò)散，因此可以看作是平面二維擴(kuò)散。重金屬元素的擴(kuò)散只能通過(guò)土壤，不考慮其他途徑（比如大氣）的。3.基本概念約定及符號(hào)說(shuō)明Igeo地質(zhì)累積指數(shù)（Muller 指數(shù)）Cij所采第i 個(gè)樣品中元素j 的濃度Bj元素j 的背景值Pi第i 個(gè)樣品的綜合污染指數(shù)Pj第j 種污染物的綜合污染指數(shù)2Iave地質(zhì)累積指數(shù)的算數(shù)平均值ImaxuDftt0地質(zhì)累積指數(shù)的最大值物質(zhì)濃度物質(zhì)的等效擴(kuò)散系數(shù) 物質(zhì)在時(shí)間內(nèi)流入時(shí)刻體積的量取樣時(shí)刻(x,y )城區(qū)內(nèi)某一為位置的坐標(biāo)（以題目中所擬規(guī)則為準(zhǔn)）h海拔高度（ h(i,j )指網(wǎng)格(i,j )處的海拔高度）Dev采樣點(diǎn)上

6、模型值與實(shí)際測(cè)量值的偏差4.問(wèn)題分析此項(xiàng)研究按照要求可以分為四個(gè)子問(wèn)題。四個(gè)子問(wèn)題逐層深入，從了解城區(qū)重金屬的空間分布狀況及不同區(qū)域污染程度入手，進(jìn)而分析出污染的主要原因，再通過(guò)分析不同污染物的在該城區(qū)的特征，從而建立污染物的模型，通過(guò)對(duì)模型的分析求解確定污染源的位置，最后，以問(wèn)題三的模型為基礎(chǔ)，收集其他信息，建立更完備的模型以研究城市地質(zhì)環(huán)境的演變模式。問(wèn)題分析導(dǎo)圖如下（圖 1）：圖 1 問(wèn)題分析導(dǎo)圖5.模型建立與求解5.1.問(wèn)題一5.1.1. 問(wèn)題分析為了解決問(wèn)題一和問(wèn)題二，需要對(duì)題目中所給數(shù)據(jù)進(jìn)行深入挖掘。由于題目中并沒(méi)有直接給出該城區(qū)的平面規(guī)劃圖，所以需要根據(jù)采樣點(diǎn)的位置坐標(biāo)和采樣點(diǎn)所

7、在的功能分區(qū)，還原城區(qū)的平面規(guī)劃圖（按功能分區(qū)），即功能區(qū)劃圖，這樣方便空間分析的進(jìn)行。然后，利用分別畫(huà)出每種元素在城區(qū)的濃度等值線圖，結(jié)合功能區(qū)劃圖可以直觀清晰的看到 8 種重金屬元素的空間分布情況。為了分析不同區(qū)域的重金屬污染程度，需要結(jié)合和該地區(qū)的實(shí)際情況建立重金屬污染程度評(píng)定模型，確定不同地區(qū)的污染級(jí)別和污染特點(diǎn)，繪制污染程度分級(jí)圖，分析后給出污染情況的結(jié)3問(wèn)題四深入研究演變模式問(wèn)題三 確定污染源問(wèn)題二分析污染原因問(wèn)題一了解污染狀況論。再結(jié)合污染物分布的等值線圖和功能區(qū)劃圖，分析出污染的主要原因。數(shù)據(jù)挖掘繪制污染物濃度等值線圖繪制功能五區(qū)圖綜合指數(shù)法地累積法污染污染程度評(píng)價(jià)污染物空間分

8、布物來(lái)源圖 2 問(wèn)題一二解題思路導(dǎo)圖5.1.2. 繪制城區(qū)功能區(qū)劃圖由于題目中所給的采樣點(diǎn)是離散的，所以，為了較準(zhǔn)確的畫(huà)出功能區(qū)劃圖，利用圖 3）：中的 griddata 函數(shù)和鄰近點(diǎn)插值（nearest）的方法，還原出城區(qū)功能區(qū)劃圖（如N圖 3 功能區(qū)劃圖由圖可見(jiàn)，山區(qū)主要分布于城區(qū)西部；工業(yè)區(qū)主要集中于東部，在城區(qū)中部也有少部分工業(yè)區(qū)；交通區(qū)分布較為廣泛，城區(qū)中心地帶較為集中；生活區(qū)零散的分布于城區(qū)4分析污染原因查閱文獻(xiàn)確定污染程度評(píng)價(jià)方法西南部和中部。5.1.3. 重金屬元素空間分布研究仍然采用等值進(jìn)行研究，即利用 griddata在研究重金屬空間分布上函數(shù)和立方插值法（cubic），繪

9、制重金屬元素在城區(qū)的等值線圖。圖 4圖 11 分別為 As、 Cd、Cr、Cu、Hg、Ni、Pb、Zn 元素的城區(qū)濃度等值線圖。圖中等值線越密集表示元素濃度越大，濃度局部峰值所在地的等值線特征為等值線閉周逐漸降低。小圓圈，濃度由中心向四圖 4 As 濃度等值線圖圖 5 Cd 濃度等值線圖圖 6 Cr 濃度等值線圖圖 7 Cu 濃度等值線圖圖 8 Hg 濃度等值線圖圖 9 Ni 濃度等值線圖圖 10 Pb 濃度等值線圖圖 11 Zn 濃度等值線圖結(jié)果與分析：5As 元素在城區(qū)的所有區(qū)域均有分布，但工業(yè)區(qū)附近濃度較高，而其他區(qū)域濃度相對(duì)較低。Cd 元素在城區(qū)的所有區(qū)域均有分布，且總體濃度較高。其中

10、，生活區(qū)、工業(yè)區(qū)和交通區(qū)中多處出現(xiàn)濃度。Cr 元素主要分布在城區(qū)的西南部，東北部的山林地區(qū)少有分布，其中在西南生活區(qū)出現(xiàn)濃度，其他分布區(qū)域濃度較低且均分布于生活區(qū)附近。Cu 元素的分布區(qū)域很少，僅占城區(qū)面積的 1/5，而且集中分布于西南生活區(qū)，其他區(qū)域分布較少，且均分布于生活區(qū)附近。Hg 元素分布區(qū)域較少，而且相互較孤立，主要集中在生活區(qū)及其周邊地區(qū)，但其濃度較高。Ni 元素在整個(gè)城區(qū)幾乎均有少量分布，濃度出現(xiàn)在了西南生活區(qū)附近，在東北山區(qū)中的交通區(qū)也出現(xiàn)了局部濃度。Pb 元素的分布區(qū)域主要集中在中部和西南部的生活區(qū)、工業(yè)區(qū)和交通區(qū)，林區(qū)分布較少，且濃度出現(xiàn)在西南部的生活、交通、工業(yè)的區(qū)域。Z

11、n 元素的分布范圍較廣，其中集中分布于中部和西南部的生活、工業(yè)、交通東北山區(qū)地帶的交通區(qū)附近也有少量分布。區(qū)，總體來(lái)看，呈現(xiàn)海拔低處污染嚴(yán)重，海拔高處污染輕；工業(yè)區(qū)，交通區(qū)污染嚴(yán)重，林地污染輕的特點(diǎn)。5.1.4. 城區(qū)重金屬污染程度分析近年來(lái)，國(guó)內(nèi)外對(duì)土壤環(huán)境價(jià)方法的研究正處于蓬勃發(fā)展階段，新的方法還在不斷出現(xiàn)，每法都有其使用范圍。在本問(wèn)題中，為了更準(zhǔn)確的對(duì)城區(qū)中不同區(qū)域的重金屬污染程度進(jìn)行分析，決定從同種元素在不同區(qū)域的污染情況（單一影響）和不同區(qū)域內(nèi)所有元算的綜合污染程度（所有綜合影響）兩個(gè)角度進(jìn)行評(píng)價(jià)。單一影響的評(píng)價(jià)，首先想到了被廣泛的單因子指數(shù)評(píng)價(jià)法，但單因子指數(shù)評(píng)價(jià)法中標(biāo)準(zhǔn)值的選擇尤

12、為重要，而本問(wèn)題是個(gè)抽象，所給數(shù)據(jù)有限，可以作為標(biāo)準(zhǔn)值的只有國(guó)家土壤環(huán)境質(zhì)量標(biāo)準(zhǔn)（GB 15618-1995），然而國(guó)標(biāo)文件中明確規(guī)定 “本標(biāo)準(zhǔn)適用于農(nóng)田、蔬菜地、茶園、果園、牧場(chǎng)、林地、自然保護(hù)區(qū)等地的土壤?！?，其標(biāo)準(zhǔn)并不適用于城區(qū)土壤評(píng)價(jià)，因此對(duì)比多種評(píng)價(jià)方法后，法。最終選擇地累積指數(shù)綜合影響的評(píng)價(jià)，考慮到了所有污染因子的平均值和方差是兩個(gè)同等重要的量值，其中任何一個(gè)值很高，都反映了土壤受到污染，只有兩個(gè)值都很低的情況下土壤才是清潔的。綜合考慮兩個(gè)量值的共同影響，想到了近些年來(lái)在環(huán)境評(píng)價(jià)方面被廣泛采用的綜合污染指數(shù)法，并對(duì)其進(jìn)行了修正。5.1.4.1. 地質(zhì)累積指數(shù)法地質(zhì)累積指數(shù)（Geoa

13、cCumulation Index）通常稱(chēng)為 Muller 指數(shù)，是 20 世紀(jì) 60 年代晚期在歐洲發(fā)展起來(lái)的廣泛用于研究沉積物及其它物質(zhì)的中重金屬污染程度的定量指標(biāo)，其表達(dá)式如下：CI logij(1)geo21.5Bj6其中，Cij 為所采第i 個(gè)樣品中元素j的修正系數(shù)。的濃度；Bj 為元素j 背景值的平均值；1.5 為背景值根據(jù) Forstner 的分級(jí)標(biāo)準(zhǔn),地質(zhì)累積指數(shù)可分為 7 個(gè)等級(jí)，0-6 級(jí)表示污染程度由無(wú)（6 級(jí)）的元素含量可能是背景值的幾百倍，不同的級(jí)別也代表著不同到極強(qiáng)，的污染程度，分級(jí)標(biāo)準(zhǔn)如表 1：表 1Muller 地質(zhì)積累指數(shù)分級(jí)標(biāo)準(zhǔn)5.1.4.2. 基于 Mul

14、ler 指數(shù)的綜合污染指數(shù)評(píng)價(jià)法對(duì)污染指數(shù)法進(jìn)行改進(jìn)，將單因子指數(shù)換成地質(zhì)累積指數(shù)，同時(shí)對(duì)其進(jìn)行加權(quán)平均，每種元素的權(quán)重按其毒性系數(shù)大小進(jìn)行確定，從而得到基于 Muller 指數(shù)的綜合指數(shù)，其數(shù)學(xué)表達(dá)式為：(W I (I )22Pj ijij max,j 1,2,(2),8i2(W I (I )22Pj ijij max,i 1,2,(3),j2 1,j 1,2,Wj(4),8.j其中，Pi 是第i 個(gè)樣品的綜合污染指數(shù)，Pj 是第j 種污染物的綜合污染指數(shù)；Iij 是第i 個(gè)樣品中元素j 的地質(zhì)累積指數(shù)；I是地質(zhì)累積指數(shù)的最大值；W 是第j 種元素maxj按期毒性系數(shù)2大小所計(jì)算的權(quán)重，權(quán)重

15、的計(jì)算方法如下：7地質(zhì)積累指數(shù) Igeo分級(jí)污染程度Igeo 00無(wú)污染0 Igeo 11輕度-中等污染1 Igeo 22中等污染2 Igeo 33中等-強(qiáng)污染3 Igeo 44強(qiáng)污染4 Igeo 55強(qiáng)-極嚴(yán)重污染5 Igeo 66極嚴(yán)重污染計(jì)算樣本的綜合污染污染指數(shù)，繪制城區(qū)重金屬元素綜合污染等值線圖（圖 12）：圖 12 城區(qū)綜合污染等值線圖由圖可見(jiàn)，整個(gè)城區(qū)重金屬污染主要集中在城區(qū)西南部，中部工業(yè)區(qū)，和南部生活區(qū)。其余區(qū)域，尤其是東北山地區(qū)域，受污染程度較輕。5.2.問(wèn)題二5.2.1. 污染物來(lái)源研究土壤重金屬污染源主要有灌溉（特別是污灌）、固體廢棄物處置（污泥、等）、和施肥以及大氣沉

16、降物等。我采與冶煉等，具體內(nèi)容見(jiàn)表 2。壤重金屬主要污染源是污灌和污泥施用、礦山開(kāi)表 2 土壤重金屬污染的主要來(lái)源8來(lái)源污染物自然環(huán)境因地而異礦產(chǎn)開(kāi)采、冶煉、加工排放的廢氣、廢水和廢渣鉻、砷、鉛、鎳煤和石油燃燒過(guò)程中排放的飄塵鉻、砷、鉛電鍍工業(yè)廢水鉻、鎘、鎳、鉛、銅電池電子工業(yè)排放的廢水、鎘、鉛工業(yè)排放的廢水化工制革工業(yè)排放的廢水鉻、鎘生活污水、工礦企業(yè)污水、混合污水、砷、鉻、鉛、鎘汽車(chē)尾氣和汽車(chē)輪胎磨損產(chǎn)生的粉塵鉛、鋅、鎘、鉻、銅固體廢棄物鎘、鉻、銅、鋅、鉛等砷、銅元素AsCdCrCuHgNiPbZn毒性系數(shù)10302540551Wij0.1020.3060.0210.0510.4080.0

17、510.0510.010圖 13 污染物功能區(qū)分布圖5.2.2. 城區(qū)重金屬污染原因分析通過(guò)對(duì)比同種元素在不同功能區(qū)平均濃度與 Mule 指數(shù)，可以確定造成該重金屬污染的主要原因（工業(yè)廢水、生活污水、汽車(chē)尾氣，固體廢棄物）。表 3 該城區(qū)不同功能區(qū)表層土壤重金屬平均濃度表 4 該城區(qū)不同功能區(qū)表層土壤重金屬 Muller 指數(shù)評(píng)價(jià)9生活區(qū)工業(yè)區(qū)山區(qū)交通區(qū)公園綠地區(qū)AsIgeo0.2160.425-0.4190.0810.213分級(jí)11011CdIgeo0.5721.011-0.3560.8850.525分級(jí)12011CrIgeo0.5700.200-0.2550.320-0.092分級(jí)1101

18、0CuIgeo1.3192.687-0.1931.6520.609分級(jí)23021HgIgeo0.8263.613-0.3583.0891.131分級(jí)14042geo-0.0090.103-0.256-0.066-0.27101000分級(jí)PbIgeo0.5721.001-0.3470.4500.385分級(jí)12011ZnIgeo1.1951.425-0.4981.2300.576分級(jí)22021As (g/g)Cd (ng/g)Cr (g/g)Cu (g/g)Hg (ng/g)Ni (g/g)Pb (g/g)Zn (g/g)生活區(qū)6.27289.9669.0249.4093.0418.3469.11

19、237.01工業(yè)區(qū)7.25393.1153.41127.54642.3619.8193.04277.93山區(qū)4.04152.3238.9617.3240.9615.4536.5673.29交通區(qū)5.71360.0158.0562.21446.8217.6263.53242.85公園綠地區(qū)6.26280.5443.6430.19114.9915.2960.71154.242015As (g/g)Cd (ng/g) Cr (g/g)10Cu (g/g)Hg (ng/g) Ni (g/g)5Pb (g/g)0生活區(qū)工業(yè)區(qū)山區(qū)交通區(qū) 公園綠地區(qū)結(jié)果分析：由表 4 分析得，在山區(qū)，8 種重金屬的 Mul

20、ler 指數(shù)均小于零，無(wú)污染，接下來(lái)只就其余四個(gè)功能區(qū)進(jìn)行比較分析。As 在生活區(qū)、工業(yè)區(qū)、交通區(qū)、公園綠地區(qū)的 Muller 指數(shù)分級(jí)均屬于 1 級(jí)，屬輕-中度污染，比較 Muller 指數(shù)：工業(yè)區(qū)生活區(qū)公園綠地區(qū)交通區(qū)，且指數(shù)值均小于 0.5，偏向于輕度污染?？梢?jiàn) As 來(lái)自于工業(yè)污染和生活污染。Cd 在工業(yè)區(qū)達(dá)到了中度污染，在生活區(qū)、交通區(qū)和公園綠地區(qū)均為輕-中度污染，Muller指數(shù)：交通區(qū)生活區(qū)公園綠地區(qū)0.5，三區(qū)均偏向于中度污染，可見(jiàn)，Cd 在四個(gè)區(qū)的污染程度相對(duì)嚴(yán)重。Cd 來(lái)自汽車(chē)尾氣污染和生活污染。Cr 在公園綠地區(qū)的 Muller 指數(shù)分級(jí)屬 0 級(jí)，無(wú)污染；在工業(yè)區(qū)、生活

21、區(qū)和交通區(qū)均屬于輕度-中度污染，比較三區(qū)的 Muller 指數(shù)：生活區(qū)0.5交通區(qū)工業(yè)區(qū)?？梢?jiàn) Cr 的污染主要來(lái)自生活廢水。Cu 在工業(yè)區(qū)的 Muller 指數(shù)分級(jí)為 3，屬中等-強(qiáng)污染，在交通區(qū)和生活區(qū)的 Muller 指數(shù)分級(jí)為 2，屬中度污染，且在交通區(qū)更嚴(yán)重，在公園綠地區(qū)則屬與輕-中度污染，總體來(lái)看，Cu 的污染強(qiáng)度較大，范圍較廣。Cu 污染來(lái)自交通區(qū)。Hg 在工業(yè)區(qū)和交通區(qū)的 Muller 指數(shù)分級(jí)均為 4 級(jí)，達(dá)到了強(qiáng)污染的程度，且在工業(yè)區(qū)更嚴(yán)重一些，在公園綠地區(qū) Muller 指數(shù)分級(jí)為 2 級(jí)，為中度污染，在生活區(qū)污染程度為輕-中等污染，綜合來(lái)看，Hg 的污染范圍集中且強(qiáng)度較

22、大。Hg 污染主要來(lái)自工業(yè)污染。Ni 的污染范圍很小，在生活區(qū)、交通區(qū)、公園綠地區(qū)均無(wú)污染，只在工業(yè)區(qū)出現(xiàn)輕-中等污染，Muller 指數(shù)僅為 0.103，偏向于輕度污染，可見(jiàn)，Ni 的污染主要來(lái)自工業(yè)區(qū)。Pb 在工業(yè)區(qū)的 Muller 指數(shù)分級(jí)為 2 級(jí)，達(dá)到了中等污染，教嚴(yán)重，在生活區(qū)、交通區(qū)和公園綠地區(qū)的污染均為輕-中度污染，比較 Muller 指數(shù)：生活區(qū)0.5交通區(qū)公園綠地區(qū)。因此Pb 污染主要來(lái)自生活污染。Zn 在生活區(qū)、工業(yè)區(qū)和交通區(qū)的 Muller 指數(shù)分級(jí)均為 2 級(jí)，已經(jīng)達(dá)到了中度污染，比較 Muller 指數(shù)：工業(yè)區(qū)交通區(qū)生活區(qū)，在公園綠地區(qū)的污染為中輕度污染（Igeo0

23、.5），可見(jiàn) Zn 的污染范圍較廣，總體程度較高。主要來(lái)自工業(yè)污染、交通污染和生活污染。5.3. 問(wèn)題三5.3.1. 問(wèn)題分析一般的濃度擴(kuò)散問(wèn)題，是根據(jù)己知源項(xiàng)和初邊值來(lái)求區(qū)域濃度場(chǎng)，此為正問(wèn)題；若源項(xiàng)的某些信息或很難直接測(cè)得，但它可以通過(guò)一些附加信息間接求得，這就了源項(xiàng)識(shí)別反問(wèn)題。問(wèn)題三要求確定污染源的具置，即屬于反問(wèn)題。圖 14 問(wèn)題三解題思路導(dǎo)圖首先建立重金屬擴(kuò)散的數(shù)學(xué)模型，由于模型中涉及到若干個(gè)參量，因本題中，此需要充分挖掘題目所給數(shù)據(jù)中的信息，對(duì)參數(shù)做一個(gè)估計(jì)；然后，在污染源可能出現(xiàn)的區(qū)域中任意選定某個(gè)位置(x0, y0 ) ，從而將問(wèn)題轉(zhuǎn)化成已知源和初值狀態(tài)的正問(wèn)題。求得到采樣點(diǎn)實(shí)

24、際測(cè)量值與模型求解值之間的偏差，讓(x0, y0 ) 遍歷所有可行解，總偏10應(yīng)用模型，找到污染源確立模型（分區(qū)域）計(jì)算模型參數(shù)（參數(shù)集總處理）建立模型（改進(jìn)的傳質(zhì)微分方程）分析特征（等值圖）差最小者即為所要求的污染源。計(jì)算最小誤差函數(shù)對(duì)源進(jìn)行搜索，從而確定污染源的位置。通過(guò)閱讀文獻(xiàn)，得到了物質(zhì)擴(kuò)散的一般擴(kuò)散方程，而確定方程的關(guān)鍵在于確定擴(kuò)散系數(shù)和阻滯因子，這兩個(gè)與所處的土壤環(huán)境、重金屬元素理化性質(zhì)、溫度、壓強(qiáng)等多種自然環(huán)境有關(guān)，通實(shí)驗(yàn)確定或由經(jīng)驗(yàn)公式確定，經(jīng)驗(yàn)公式適用范圍有限卻條件苛刻，而本題是一個(gè)抽象模型，且研究對(duì)象多樣，難以找到適用的經(jīng)驗(yàn)公式，當(dāng)然，更不具備實(shí)驗(yàn)測(cè)定參數(shù)的條件。所以，決定

25、利用反分析方法確定模型參數(shù)。由于擴(kuò)散系數(shù)和阻滯因子存在著確定的關(guān)系，根據(jù)工程師原則，利用集數(shù)的，考慮在將擴(kuò)散系數(shù)和阻滯因子兩個(gè)參數(shù)化作一個(gè)參數(shù)，即等效擴(kuò)散系數(shù)，這樣可以較為方便的進(jìn)行計(jì)算求解。5.3.2. 模型 I首先，建立一個(gè)簡(jiǎn)化的模型，不考慮地勢(shì)對(duì)元素?cái)U(kuò)散造成的影響根據(jù)傳質(zhì)學(xué)基本原理，類(lèi)比熱傳導(dǎo)方程，建立重金屬污染在土壤中擴(kuò)散的微分方程如下：u D2u f(5)t其中，2 為算符；u 為元素濃度，且為 x, y,t 的函數(shù)；D 為等效擴(kuò)散系數(shù)，f 為時(shí)間外界重金屬流入量。將式(5)寫(xiě)成分量的形式則為：u(x, y,t) 2u(x, y,t)2u(x, y,t)f (x, y t)D(6)t

26、x2y2在本問(wèn)題中，首先將污染源看作一個(gè)點(diǎn)，不妨設(shè)污染源所在點(diǎn)坐標(biāo)為(x ,y )。00那么，f (x, y, t) f0 ,(x, y) (x0 , y0 )(7)0, (x, y) (x , y )00f (x, y, 0) BE其中， BE 為土壤中重金屬元素的背景值。這是一個(gè)普適方程，實(shí)際上，D 會(huì)隨環(huán)境、物質(zhì)濃度等(8)的變化而發(fā)生變化，而基于對(duì)問(wèn)題一的分析已經(jīng)確定了污染源所在的大致區(qū)域，在局部范圍內(nèi)，環(huán)境等因素基本相同，所以在局部進(jìn)行問(wèn)題研究可以將 D 取作常數(shù)。5.3.3. 模型 II為了便于計(jì)算機(jī)進(jìn)行求解，模型一中的連續(xù)方程進(jìn)行離散化處理，建立差分方程模型。以步長(zhǎng)s 將x 和y

27、 劃分為均勻的網(wǎng)格，任意網(wǎng)格(i,j )，可以與其正對(duì)的 4個(gè)網(wǎng)格之間進(jìn)行物質(zhì)傳遞，在此模型中仍然假定擴(kuò)散系數(shù) D 是個(gè)與 x, y 取值無(wú)關(guān)的11常數(shù)。用三點(diǎn)中心二階差分來(lái)近似u 的二階導(dǎo)數(shù)，有：2u u(xi1, y j ,tk ) 2u(xi , y j ,tk ) u(xi1, y j , tk )(9)x2則差分方程如下：s2) D(u(xi1, y j ,tk ) 2u(xi , y j ,tk ) u(xi1, y j ,tk )u(x , y ,t) u(x , y ,tij k 1ij ks2u(xi , y j1,tk ) 2u(xi , y j , tk ) u(xi ,

28、 y j 1, tk ) f (x , y,t )(10)ij ks2其中，f (x , y , t) f0 ,當(dāng)(x0 , y0 )屬于網(wǎng)格(i, j)ij0, otherwisek(11) f (x , y , 0) BEij5.3.4. 模型 III模型 III 引入海拔的影響。考慮網(wǎng)格(i,j )與(i 1,j )之間的物質(zhì)擴(kuò)散，兩者的濃度差為u(xi1, yj ) u(xi , yj ) ，濃度差為正，表示物質(zhì)從(i 1, j) 流向(i, j) ，(i, j) 濃度增加，否則反之；若兩者的相對(duì)海拔差異較大，則會(huì)對(duì)兩者之間的物質(zhì)擴(kuò)散起促進(jìn)或者阻礙作用：表 5 海拔對(duì)物質(zhì)擴(kuò)散的影響由上

29、表可知，記相對(duì)高度對(duì)擴(kuò)散速度的促進(jìn)或阻礙作用為 (i, j,i 1, j) ，則12濃度比較海拔高度比較高度對(duì)擴(kuò)散速度影響(i, j) 的濃度變化u(i, j) u(i 1, j)h(i, j) h(i 1, j)促進(jìn)流入h(i 1, j) (u(i 1, j) u(i, j) h(i, j)h(i, j) h(i 1, j)阻礙流入h(i 1, j) (u(i 1, j) u(i, j) h(i, j)u(i, j) u(i 1, j)h(i, j) h(i 1, j)促進(jìn)流出h(i, j) (u(i 1, j) u(i, j) h(i 1, j)h(i, j) h(i 1, j)阻礙流出h

30、(i, j) (u(i 1, j) u(i, j) h(i 1, j) h(i 1,j) , u(i, j) u(i 1, j)h(i, j)h(i, j) (i, j, i 1, j) (12), u(i, j) u(i 1, j) h(i 1, j)擴(kuò)散方程可改寫(xiě)為：u(xi , y j , tk 1 ) u(xi , y j , tk ) D (i, j, i 1, j) u(x, y , t ) u(x , y , t ) (i, j, i 1, j) u(x, y , t ) u(x , y , t )i1ji 1jkij kkij ks2(13) (i, j, i, j 1) u(

31、x , y, t ) u(x , y , t ) (i, j, i, j 1) u(x , y, t ) u(x , y , t )j 1 kj 1 kiij kiij kf (xi , y j , tk )f (x , y , t) f0 ,當(dāng)(x0 , y0 )屬于網(wǎng)格(i, j)0, otherwiseij k f (x , y , 0) BEij5.3.5. 模型求解5.3.5.1. 離散化過(guò)程從等值線圖中可以看出，一個(gè)污染源影響的范圍是有限的，在一定距離外，可以視作不受影響；可以用一個(gè)矩形限定影響的范圍，從而減小解空間。所以從污染源產(chǎn)生的污染，可以看作完全流入到該矩形區(qū)域，即f t

32、d u(x, y)dxdy(14)5.3.5.2. 參數(shù)估計(jì)取樣時(shí)間t本題中，記初始時(shí)刻t 0 ，樣本數(shù)據(jù)的時(shí)刻t0 未知，且與時(shí)間有關(guān)；當(dāng)變換時(shí)間變量，如令t ，則方程(5)變?yōu)閡 2u2uD() f(15)x2y2僅是系數(shù)同比例變化。因此，不妨認(rèn)定取樣時(shí)刻t 1，即題中所給的測(cè)量值是t 1時(shí)的數(shù)據(jù)。時(shí)間重金屬流入量f因?yàn)槿 1 ，所以， f 理論值計(jì)算公式為：f d u(x, y)dxdy(16)13其中， 1340g / m3 為土壤密度； d 為土壤厚度，題目給出取表層土壤 0.01m。由于u(x, y) 未知，只知道幾個(gè)離散采樣點(diǎn)上濃度值，通過(guò)觀察數(shù)據(jù)發(fā)現(xiàn)，所有的污染區(qū)域，都是由若

33、干個(gè)污染程度較低、接近于背景值的采樣點(diǎn)和一個(gè)或兩個(gè)靠近污染區(qū)域幾何中心、污染程度嚴(yán)重的采樣點(diǎn)所組成，因此用 f d S umax 來(lái)估計(jì) f ，其中 S 為3污染區(qū)域的面積， umax 為靠近中心采樣點(diǎn)的重金屬含量值。直觀上來(lái)講，就是用棱錐的體積來(lái)估計(jì)曲面u(x, y) 包圍的體積（如圖 15）。圖 15 曲面體積估算示意圖雖然理論上可以令t 1 ，但是由于計(jì)算機(jī)模擬擴(kuò)散與實(shí)際擴(kuò)散并不完全符合，以及參數(shù)D 需要在小范圍內(nèi)搜索確定，t 的取值也會(huì)影響到最終解得的污染源位置。當(dāng)總輸入量恒定時(shí)( f(t )dt 為常數(shù))， f 可以是：輸入，衰減輸入等。輸入：實(shí)例的計(jì)算機(jī)仿真均采用該種輸入法，其特點(diǎn)

34、是簡(jiǎn)單，計(jì)算速度快，擴(kuò)散方式呈現(xiàn)線性擴(kuò)散，即從污染源開(kāi)始，污染程度線性減少。輸入（t 1 ），持續(xù)穩(wěn)定持續(xù)穩(wěn)定輸入：相比輸入，該輸入方式更符合生活生產(chǎn)實(shí)際，但需要多確定一個(gè)參數(shù)t ，即持續(xù)輸入時(shí)間。在t 確定和污染物輸入總量確定的情況下，也比較簡(jiǎn)單，計(jì)算速度比較快。衰減輸入：在特定情況下，污染源可能慢慢減少污染量，如溢油封堵。需要給出f f(t )dt =污染物總量，計(jì)算較為復(fù)雜，速度也較慢。與時(shí)間的函數(shù)，可以根據(jù)不同的污染源特點(diǎn)，選擇不同的污染輸入方式，使得模型求解結(jié)果更加精確。擴(kuò)散系數(shù)D14擴(kuò)散系數(shù)D 只與土壤本身性質(zhì)有關(guān)。由于樣本點(diǎn)過(guò)少，所給信息，難以直接從所給采樣數(shù)據(jù)中得到；同時(shí)經(jīng)過(guò)試

35、探性計(jì)算發(fā)現(xiàn)D 的取值變化不大，在 0.2 附近，所以對(duì)D 的估計(jì)可以在求解污染源的過(guò)程中，讓D 從 0,1 枚舉到 0.3，與實(shí)際情況最為匹配的值視為對(duì)D 的估計(jì)D。5.3.5.3. 污染源確定若已知參數(shù)D ,f ，那么，在污染區(qū)域取任意點(diǎn)x0 ,y0 作為污染源，方程(5)便成為已知初始狀態(tài)和各個(gè)參變量的擴(kuò)散正問(wèn)題，進(jìn)而可以求出采樣時(shí)刻（即t 1 時(shí)刻）第i個(gè)采樣點(diǎn)位置處的濃度值u(x(i), y(i),1) ，與實(shí)際測(cè)量值u0 (i) 之間比較，可用偏差kDev (u(x(i), y(i),1) u (i)2,i 1, 2,., k(17)0i1來(lái)衡量。使得 Dev 最小的(x0 , y

36、0 ) 即為題目所求污染源。5.3.5.4. 算法分析f 進(jìn)行估計(jì)，估計(jì)值為 f 。ax ymin, ymax ；對(duì)第一步：確定污染范圍xm第二步：任取污染范圍中的一點(diǎn)(x0, y0 ) 作為污染源，令 D 從 0.1 變到 0.3，對(duì)于每一個(gè) D ，求解該已知源和初值狀態(tài)的正問(wèn)題，計(jì)算采樣點(diǎn)上理論值和實(shí)際值偏差 Dev ，找到使得 Dev 最小的 D ，作為 D 的估計(jì)。第三步：讓(x0 , y0 ) 遍歷所有可能的污染范圍中的點(diǎn)，比較不同(x0 , y0 ) 的 Dev ,找到使得 Dev 最小的(x0 , y0 ) ，即為所求污染源位置。5.3.5.5. 結(jié)果分析首先，根據(jù)問(wèn)題一中用插值

37、法繪制的元素分布等值線圖將污染區(qū)進(jìn)行劃分，再結(jié)合污染區(qū)的局部采樣數(shù)據(jù)，利用計(jì)算機(jī)進(jìn)行污染源搜索，從而確定每個(gè)污染區(qū)內(nèi)污染源的準(zhǔn)確位置。下面，以 Cr 元素為例，對(duì)計(jì)算機(jī)求解的結(jié)果進(jìn)行分析：15圖 16 Cr 元素污染分區(qū)圖表 6 污染源位置坐標(biāo)在圖 16 上Cr 元素污染較嚴(yán)重的三個(gè)區(qū)域分別標(biāo)記為， A(Cr ) 、B(Cr ) 、C(Cr ) 。從圖中可以判斷污染源的大致區(qū)域，但是無(wú)法準(zhǔn)確判斷，尤其像C這樣中心面積(Cr )較大的區(qū)域更加難以確定。而通過(guò)對(duì)上述模型的求解，可以較準(zhǔn)確的得到污染源所在位置的坐標(biāo)，其中， A 區(qū)域的污染源位于（21910,11890）位置附近， B 區(qū)域的污(Cr

38、 )(Cr )染源位于（10500，5930）位置附近，C(Cr ) 區(qū)域的污染源位于（3630，5808）位置附近。其他元素污染源確定結(jié)果的分析方法與 Cr 類(lèi)似，僅把結(jié)果如下（其中：左側(cè)為相應(yīng)元素污染區(qū)劃圖， x 為污染源橫坐標(biāo)值， y 為污染源縱坐標(biāo)值）：0016污染區(qū)A(Cr )B(Cr )C(Cr )x021910105003630y01189059305808表 7 As 元素污染源確定結(jié)果表表 8 Cd 元素污染源確定結(jié)果表表 9 Cu 元素污染源結(jié)果確定表表 10 Hg 元素污染源結(jié)果確定表17污染區(qū)A(Hg )B(Hg )C(Hg )x015560138003110y0920

39、026502600污染區(qū)A(Cu )B(Cu )x035402560y058003390污染區(qū)A(Cd )B(Cd )C(Cd )D(Cd )E(Cd )F(Cd )x021060176606570363040102280y011320392010050574026502960污染區(qū)A(As )B(As )C(As )D(As )E(As )x01810012530703047901320y010440,2600776077302480表 11 Ni 元素污染源結(jié)果確定表表 12 Pb 元素污染源結(jié)果確定表表 13 Zn 元素污染源結(jié)果確定表5.3.5.6. 模型檢驗(yàn)（相關(guān)性分析）對(duì)比分析污染

40、區(qū)劃圖，猜想A和A，C、B和B，F(xiàn)和C， A(Cr )(Ni )(Cr )(Ni )(Pb )(Cd )(Pb )(Cu )和E Zn ， B(Cu ) 和G(Zn ), E(As ) 和F(Cd ) 可能有相同的污染源。用 spss相應(yīng)重金屬元素進(jìn)行相關(guān)性分析，分析結(jié)果如下：分別對(duì)這幾個(gè)區(qū)域的18污染區(qū)A(Zn )B(Zn )C(Zn )D(Zn )E(Zn )F(Zn )x014160128409520512036204250y0947034704580688057302680污染區(qū)A(Pb )B(Pb )C(Pb )D(Pb )x014140437014703910y0123704968

41、23002590污染區(qū)A(Ni )B(Ni )x0219303540y0120205890表 14 不同元素不同污染區(qū)污染物擴(kuò)散狀況相關(guān)性分析表根據(jù) spss的相關(guān)性分析，Cr 與 Ni 在 A(Cr)區(qū)域的濃度分布高度相關(guān)，可以說(shuō)明在該區(qū)域，Cr 與 Ni 的污染源相同；基于模型計(jì)算的地理坐標(biāo)分別為（21910,11890）和（21930,12020），兩者在誤差范圍內(nèi)基本重合。另外三組分析方法同上，且模型計(jì)算的地理坐標(biāo)在誤差范圍內(nèi)基本重合。相關(guān)性分析結(jié)果對(duì)模型進(jìn)行了檢驗(yàn)，表明模型是合理的。6.模型評(píng)價(jià)與改進(jìn)6.1.模型評(píng)價(jià)6.1.1. 模型優(yōu)點(diǎn)題目要求利用取樣點(diǎn)數(shù)據(jù)，求解擴(kuò)散反問(wèn)題，確定污

42、染源位置。本題最突出的特點(diǎn)是所提供的信息非常有限，因此建立模型時(shí)，立足點(diǎn)就是如何盡可能的利用有限的數(shù)。根據(jù)模型 III 求解據(jù)，充分挖掘出數(shù)據(jù)中的隱含信息，求解出一個(gè)盡可能精確的的結(jié)果表明，所確定污染源的精度還是令人滿(mǎn)意的，而且對(duì)部分地區(qū)元素分布進(jìn)行相關(guān)性檢驗(yàn)，也在一定程度上對(duì)模型進(jìn)行了驗(yàn)證。因此總的來(lái)說(shuō)，模型 III 是合理的，也是滿(mǎn)足題目需要的。6.1.2. 模型缺點(diǎn)模型 III 的缺陷也是顯而易見(jiàn)的。土壤中重金屬的擴(kuò)散本身是個(gè)極其復(fù)雜的過(guò)程，有多個(gè)參量，涉及到擴(kuò)散、對(duì)流、吸附、降解等諸多物理過(guò)程，需要流體力學(xué)、材料學(xué)、土壤學(xué)、環(huán)境學(xué)等多個(gè)領(lǐng)域相關(guān)專(zhuān)業(yè)知識(shí)。但是受限于數(shù)據(jù)信息的匱乏，許多影

43、響只能夠化簡(jiǎn)，合并，甚至忽略不計(jì)。比如未考慮重金屬在土壤中的降解，模型 III 中相對(duì)海拔對(duì)擴(kuò)散速率的影響因子 (x, y) ，實(shí)際上是對(duì)水生的對(duì)流影響的粗略估計(jì)。又如參數(shù) D ，和土壤參數(shù)有關(guān)，可以用實(shí)驗(yàn)的方法求出精確值；但在本題中，卻需要用其他數(shù)據(jù)來(lái)估計(jì) D ，從而會(huì)累計(jì)誤差。如果能有較為充分的數(shù)據(jù)，那么可以在模型 III的基礎(chǔ)上，建立起更符合土壤中重金屬的擴(kuò)散模式、更符合實(shí)際情況的模型。6.2.信息收集信息的收集需要考慮信息的有效性和可獲取性。有效性是指信息有助取得更好的結(jié)果，可獲取性是指信息是易于獲得的。例如，模型中，若能得到污染源時(shí)間內(nèi)重金與t 的比值，進(jìn)而直接得到屬的排放量（即可以

44、得到模型 III 中ff），必然可以使得結(jié)果更加精確；但在實(shí)際情況，首先污染源未知，排放量的統(tǒng)計(jì)更是非常，甚至還會(huì)出現(xiàn)相關(guān)故意的現(xiàn)象，因此污染源時(shí)間內(nèi)重金屬的排放量這個(gè)信息是不易19污染區(qū)A(Cr)-A（Ni）C(Gr)-B(Ni)F(Cd)-C(Pb)A(Cu)-E(Zn)Pearson系數(shù)0.9980.8100.9880.892相關(guān)程度高度相關(guān)高度相關(guān)高度相關(guān)高度相關(guān)獲取的?？紤]到以上兩點(diǎn)，再結(jié)合實(shí)際中重金屬的擴(kuò)散方式，可以通過(guò)采樣、實(shí)驗(yàn)獲取以下信息：Dk擴(kuò)散引起的重金屬擴(kuò)散系數(shù)，可由實(shí)驗(yàn)得出；對(duì)流引起的重金屬擴(kuò)散系數(shù)，可由實(shí)驗(yàn)得出；重金屬衰減系數(shù)，可由實(shí)驗(yàn)得出；v(x, y,t)在點(diǎn)(x

45、, y) 出的水流速度和方向，可通過(guò)實(shí)地測(cè)出；6.3.模型改進(jìn)土壤污染在土壤中隨著土壤水分的運(yùn)移而運(yùn)動(dòng)，同時(shí)土壤溶質(zhì)在運(yùn)移的過(guò)程中還會(huì)受到土壤基質(zhì)、土壤溶質(zhì)本身的物理化學(xué)性質(zhì)以及土壤中其他化學(xué)成分的影響。土壤重金屬濃度發(fā)生變化的物理機(jī)理主要包括對(duì)流，擴(kuò)散以及因土壤自?xún)裟芰椭参镂辗纸庠斐傻乃p。考慮x y 的小面元在t 時(shí)間里的重金屬質(zhì)量m 變化，擴(kuò)散作用：擴(kuò)散是由重金屬原子的熱運(yùn)動(dòng)引起的，具有從高濃度處向低濃度處運(yùn)移，以求達(dá)到濃度相等的趨勢(shì)。與熱擴(kuò)散方程類(lèi)似，由 Fourier 實(shí)驗(yàn)定理，擴(kuò)散運(yùn)動(dòng)引起的小面元內(nèi)，2u2u重金屬質(zhì)量變化量為m1 D( x2 y2 )xyt對(duì)流作用對(duì)流是指土壤

46、重金屬隨著土壤水分的運(yùn)移過(guò)程而進(jìn)行的運(yùn)動(dòng)過(guò)程。，水流方向 v 與重金屬濃度梯度方向n (u , u ) 之間存在角度，因此對(duì)流作用引起的小面元內(nèi)x y重金屬質(zhì)量變化量為m v (u , u )xyt (vu vu )xyt ，其中v , v 分別2x yxxyyxy是 v 的 x, y 軸分量。衰減考慮小面元內(nèi)重金屬因土壤自?xún)裟芰Χ鴾p少的質(zhì)量，一般而言，質(zhì)量衰減的速率與土壤濃度成正比，因此m3 kuxyt利用m m1 m2 m3 ，再加上外界重金屬的流入，從而2u2uuum (D( x2 y2 ) (vx x vy y ) ku f )xyt20u 2u2uuu化為微分形式即為：D() (v

47、v) ku ftx2y2x xy y與模型 III 的處理類(lèi)似，將之轉(zhuǎn)化為差分方程便可求解。參考文獻(xiàn)1等，徐州市九里礦區(qū)土壤重金屬插值分析及污染評(píng)價(jià)，測(cè)繪科學(xué)，35（5）：166-169 頁(yè),2010.231（2）， 341989.啟等，潛在生態(tài)危害指數(shù)法評(píng)價(jià)中重金屬毒性系數(shù)計(jì)算，環(huán)境科學(xué)與技術(shù)，112-115,2008.，2000.：科學(xué)，大學(xué)生數(shù)學(xué)建模競(jìng)賽輔導(dǎo)，湖南：湖南教育，環(huán)境中微量重金屬元素的污染危害與遷移轉(zhuǎn)化，5 Cleve.B.Moler，2010數(shù)值計(jì)算，：機(jī)械工業(yè)62008.7士，土壤中重金屬污染源研究進(jìn)展，有色金屬Vol.60,No.4，重金屬污染物在土壤環(huán)境系統(tǒng)中運(yùn)移的建

48、模與仿真，廣西師范大學(xué)碩,2006.8，城市不同功能區(qū)土壤重金屬分布初探，土壤學(xué)報(bào)，Vol.42，吳新民，No.3，2005.9 張長(zhǎng)波，等，污染場(chǎng)地土壤重金屬含量的空間變異特征及其污染源識(shí)別指示意義，土壤，Vol.38，No.5，2006.21附錄一：/* 污染源位置檢測(cè) */* 2011/09/11*/#include #include #include#include using namespatd;double u200200; / 污染程度矩陣double h200200; / 海拔矩陣doub污染背景值/ 樣本點(diǎn)坐標(biāo)/ 樣本點(diǎn)污染程度/ 矩陣行數(shù),列數(shù)y,樣本點(diǎn)數(shù)量kn202;do

49、uble knu10; n,m,l;double mini,minj,mind,minn=1e9;/最優(yōu)解double d; double f200;bool flag200200;/ 等效擴(kuò)散系數(shù)/ 污染物輸入量標(biāo)記bool inq200200; / 隊(duì)列內(nèi)標(biāo)記double minu200200; / 最優(yōu)解下擴(kuò)散結(jié)果矩陣double t;/ 時(shí)間/ 檢查解矩陣是否為全正bool check_itive(double u200200)for (i=0; in; i+)for (j=0; jm; j+)if (uij 0) return false;return true;/ 擴(kuò)散模擬void calc(double s;x,y)/ 污染擴(kuò)散量vectorqx,qy;/ 擴(kuò)散隊(duì)列for (i=0; in; i+) / 污染程度預(yù)設(shè)為背景值for (j=0; jm; j+)uij = back;i=0; it; i+)for (22uxy=uxy+fi; / 按照時(shí)間依次