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

1、城市表層土壤重金屬污染分析摘 要隨著全球經(jīng)濟(jì)化的迅速發(fā)展，人類活動(dòng)對城市地質(zhì)環(huán)境影響日益顯著，重金屬污染物通過各種途徑進(jìn)入土壤，土壤重金屬污染造成了生態(tài)系統(tǒng)嚴(yán)重破壞。模型一：將采樣點(diǎn)的空間坐標(biāo)與重金屬污染物濃度綜合起來，繪制出元素的空間分布圖（三維圖和等高圖），綜合使用空間分布圖，可直觀的看出各功能區(qū)的密集程度以及各種重金屬污染物濃度。通過地質(zhì)累積指數(shù)、內(nèi)梅羅指數(shù)可求出每種元素在每個(gè)區(qū)域內(nèi)的綜合污染指數(shù)，利用各種重金屬元素的綜合污染指數(shù)，建立重金屬污染物對環(huán)境污染程度的模糊綜合評判矩陣，并運(yùn)用貼近度算法，求出各功能區(qū)的相對污染程度，得到各區(qū)域污染程度大小序列：交通區(qū)工業(yè)區(qū)生活區(qū)公園綠地區(qū)山區(qū)。

2、模型二：利用因子分析法和主成分分析法，求出各金屬元素的相關(guān)系數(shù)矩陣，求得兩個(gè)主成分：對第一個(gè)主成分貢獻(xiàn)較大的元素有hg、cu，對第二主成分貢獻(xiàn)較大的元素有as 、cr 、ni。最終得出結(jié)論：該城區(qū)的污染主要是由人為活動(dòng)影響的。交通污染和工業(yè)污染是主要原因。 根據(jù)各功能區(qū)各種重金屬元素的綜合污染指數(shù)，繪制出重金屬元素污染指數(shù)直方圖，可直觀地觀察出在不同區(qū)域內(nèi)其主要污染元素，并尋找其源頭及造成因素。 模型三：綜合考慮重金屬污染物的傳播特征，根據(jù)重金屬污染的不同傳播特征，建立兩個(gè)不同的污染源確定模型。一、根據(jù)污染物傳播特征（高濃度向低濃度擴(kuò)散），判斷污染源可能在濃度較高的位置。利用各重金屬元素的污染

3、指數(shù)生成歸一化矩陣，按照各污染物的綜合污染指數(shù)，將污染區(qū)域的位置分成三個(gè)模塊。利用fminsearch函數(shù)最優(yōu)解算法和初始模糊猜測值求目標(biāo)函數(shù)最優(yōu)解對應(yīng)的空間坐標(biāo)。從而得到模糊最優(yōu)污染源位置。二、利用城區(qū)地貌特征和污染物傳播特征（高海拔向低海拔擴(kuò)散），判斷原始污染源可能在海拔較高的位置，運(yùn)用曲線擬合和微分方程的性質(zhì)，可得到多個(gè)近似污染源；對重金屬污染源位置分類討論，得出原始污染源：（15600，15286，118）。模型四：通過對模型三的綜合評價(jià)，可看出該模型缺乏重金屬污染源的分布隨時(shí)間而產(chǎn)生的變化，及其對污染源位置確定的影響。采取按一定時(shí)間梯度在整個(gè)城區(qū)進(jìn)行取樣考察的方法，對不同時(shí)間城市地貌

4、進(jìn)行插值逼近，求出城市地質(zhì)環(huán)境中重金屬污染物以及污染源在不同時(shí)間內(nèi)的分布規(guī)律，進(jìn)一步的探索人類活動(dòng)影響下城市地質(zhì)環(huán)境的演變模式。關(guān)鍵詞：內(nèi)梅羅指數(shù)、模糊綜合評判、主成分分析一、問題的重述隨著城市經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷增加，人類活動(dòng)對城市環(huán)境質(zhì)量的影響日顯突出。對城市土壤地質(zhì)環(huán)境異常的查證，以及如何應(yīng)用查證獲得的海量數(shù)據(jù)資料開展城市環(huán)境質(zhì)量評價(jià)，研究人類活動(dòng)影響下城市地質(zhì)環(huán)境的演變模式，日益成為人們關(guān)注的焦點(diǎn)。按照功能劃分，城區(qū)一般可分為生活區(qū)、工業(yè)區(qū)、山區(qū)、主干道路區(qū)及公園綠地區(qū)等，分別記為1類區(qū)、2類區(qū)、5類區(qū)，不同的區(qū)域環(huán)境受人類活動(dòng)影響的程度不同?，F(xiàn)對某城市城區(qū)土壤地質(zhì)環(huán)境進(jìn)行調(diào)

5、查。為此，將所考察的城區(qū)劃分為間距1公里左右的網(wǎng)格子區(qū)域，按照每平方公里1個(gè)采樣點(diǎn)對表層土（010 厘米深度）進(jìn)行取樣、編號，并用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種主要重金屬元素的背景值。下面我們通過建立數(shù)學(xué)模型解決以下問題：(1) 給出8種主要重金屬元素在該城區(qū)的空間分布，并分析該城區(qū)內(nèi)不同區(qū)域重金屬的污染程度。(2)

6、 通過數(shù)據(jù)分析，說明重金屬污染的內(nèi)在因素和人為因素。(3) 分析重金屬污染物的傳播特征，由此建立模型，確定污染源的位置。(4) 分析問題（3）所建立模型的優(yōu)缺點(diǎn)，為更好地研究城市地質(zhì)環(huán)境的演變模式，還應(yīng)收集哪方面信息，以及如何利用這些信息，建立更優(yōu)化的模型來確定污染源的位置。二、問題的分析由樣本點(diǎn)的坐標(biāo)，各種重金屬元素的濃度，可以繪制出三維圖形。通過內(nèi)梅羅指數(shù)方程，確定各濃度污染物的權(quán)重，結(jié)合數(shù)據(jù)可確定各重金屬污染物對環(huán)境影響的綜合污染指數(shù)，并對內(nèi)梅羅綜合污染指數(shù)貼近度模糊評判，從而確定不同區(qū)域的重金屬污染物對環(huán)境的污染程度。對各種重金屬元素進(jìn)行因子分析和主成分分析，確定各城區(qū)重金屬污染主要由

7、某幾種元素造成，具體分析該元素濃度大的原因。觀察地貌圖結(jié)合重金屬污染物的傳播特征，將其分為兩方面進(jìn)行考慮（即高濃度向低濃度擴(kuò)散：求濃度較高點(diǎn)為污染源；高海拔向低海拔擴(kuò)散：求海拔較高點(diǎn)）?？梢詫?shù)據(jù)進(jìn)行擬合，或者運(yùn)用fminsearch函數(shù)求最優(yōu)解求出污染源的大概最優(yōu)位置。研究城市地質(zhì)環(huán)境的演變模式，需要結(jié)合時(shí)間問題來綜合判斷污染源的位置。一定程度上反映出人類活動(dòng)對城市地質(zhì)環(huán)境演變的影響。三、模型的假設(shè)與符號說明3.1、 模型的假設(shè)（1）假設(shè)一：測量數(shù)據(jù)準(zhǔn)確，不考慮測量儀器的誤差及人為因素；（2）假設(shè)二: 采集數(shù)據(jù)在同一時(shí)間點(diǎn)采集；（3）假設(shè)三：地質(zhì)結(jié)構(gòu)穩(wěn)定，在短時(shí)間內(nèi)不發(fā)生改變；（4）假設(shè)四：

8、各種重金屬元素的傳播特征相同；（5）假設(shè)五：忽略生物因素對重金屬元素傳播造成的影響；3.2、 符號說明第i種元素在j功能區(qū)的地質(zhì)累積指；第i種元素在j功能區(qū)的濃度；第i種元素的背景濃度；第i種元素在第j功能區(qū)的綜合污染指數(shù)；地質(zhì)累積指數(shù)算術(shù)平均值； 地質(zhì)累積指數(shù)最大值; a歸一化后各元素污染指數(shù)矩陣；y歸一化后同樣本點(diǎn)各元素污染指數(shù)和；x分別表示x,y,海拔高度；四、模型的建立與求解4.1.問題（1）:為了描繪出8種主要重金素元素在該城區(qū)的空間分布，分析了每個(gè)采樣點(diǎn)各種重金屬元素濃度，并利用多元統(tǒng)計(jì)分析方法研究不同功能區(qū)與污染物濃度之間的相互聯(lián)系,從而得出空間累積規(guī)律。將附件1與附件2整合到一

9、起，為shuju表。對于該問題我們采用了畫不同元素在整個(gè)區(qū)域的空間分布圖，將采取樣點(diǎn)的x,y作為x軸y軸坐標(biāo)，z軸為重金屬元素在各個(gè)樣點(diǎn)的濃度大小，繪制三維立體圖；由于每種元素對于五個(gè)功能區(qū)的影響因素不同，我們可以根據(jù)shuju表，將五個(gè)區(qū)用不同的符號表示在圖中，由每種元素在各功能區(qū)的濃度及地理位置，可以畫出一個(gè)二維圖。利用matlab軟件（程序代碼見附錄表1），繪制出8種主要重金屬元素在該城區(qū)的空間分布圖，如下：根據(jù)matlab軟件畫出的重金屬元素在該城區(qū)內(nèi)的空間分布圖，可以分析出在某種重金屬元素高對應(yīng)的地理位置功能區(qū)分類的密集程度，判斷出單元素對某個(gè)功能區(qū)的污染最大。我們將數(shù)據(jù)定量化分析，

10、建立重金屬對環(huán)境的模糊綜合評判模型。 模型一：重金屬對環(huán)境污染的模糊綜合評判模型每種元素在每個(gè)功能區(qū)的綜合污染指數(shù)：首先，我們?yōu)榱搜芯砍练e物及其他物質(zhì)中重金屬污染程度的定量指標(biāo)，將shuju表按功能區(qū)序號進(jìn)行排序，令aj為第j區(qū)樣點(diǎn)各元素的濃度數(shù)據(jù)表，利用重金屬污染物的地質(zhì)累積指數(shù)(geoaccumulationindex又稱為muller指數(shù))，其表達(dá)公式如下： 其中，第i種元素在j功能區(qū)的地質(zhì)累積指數(shù)；第i種元素在j功能區(qū)的濃度；元素i的背景濃度；計(jì)算出as、cd、cr、cu、hg、ni、pb、zn八種重金屬元素分別對每個(gè)功能區(qū)的地質(zhì)積累指數(shù)。其次，突出環(huán)境要素中濃度最大的污染物對環(huán)境質(zhì)量

11、的影響，采用地質(zhì)累積指數(shù)與內(nèi)梅羅指數(shù)相結(jié)合的綜合指數(shù)法進(jìn)行污染評價(jià)，得到新的綜合指數(shù)，公式如下：其中，第i種元素在第j功能區(qū)的綜合污染指數(shù)；地質(zhì)累積指數(shù)算術(shù)平均值； 地質(zhì)累積指數(shù)最大值;matlab運(yùn)行（程序代碼見附錄2），求得每種元素在每個(gè)功能區(qū)的綜合污染指數(shù)cj,元素指數(shù)as cdcr cu hg nipbzn c10.7721.72872.8312.65552.40050.589552.37123.4138c21.43481.84661.85175.04235.73740.832162.32982.8744c30.816940.839191.37881.30471.46121.45390

12、.968550.91522c41.75422.20053.0464.38475.84392.08951.40343.6941c50.792991.69770.751682.03843.30720.523231.62172.6496采用貼近度的方法對五個(gè)功能區(qū)進(jìn)行污染程度排序：以絕對污染程度（由同一元素五個(gè)功能區(qū)的污染程度的最大值作為該元素對環(huán)境的污染程度）和五個(gè)不同地區(qū)的不同元素的污染程度作為指標(biāo)，建立模糊評價(jià)矩陣:采用貼近度來計(jì)算五個(gè)功能區(qū)的相對污染程度：生活區(qū) 工業(yè)區(qū) 山區(qū) 交通區(qū) 公園綠地區(qū) 0.6115 0.8515 0.3601 0.9721 0.4789 所得數(shù)據(jù)恰好與元素空間分布

13、圖的分布狀況相吻合，對八種元素在五個(gè)功能區(qū)的綜合污染從強(qiáng)到弱排序：交通區(qū)、工業(yè)區(qū)、生活區(qū)、公園綠地區(qū)、山區(qū)。4.2.問題（2）重金屬的主要原因：重金屬以各種化學(xué)狀態(tài)或化學(xué)形態(tài)存在于自然界中，隨著大氣、水、土壤、生物體的流動(dòng)進(jìn)入生態(tài)系統(tǒng)循環(huán)，然而重金屬循環(huán)中不能被生物降解，可在土壤中富集，所以生活中的重金屬含量會(huì)越來越多。重金屬過量存在嚴(yán)重的損壞了人們的身體健康，因此分析某種重金屬污染的主要原因十分重要。本文對重金屬污染的主要原因從兩個(gè)標(biāo)準(zhǔn)來進(jìn)行評判。首先，對于整個(gè)城區(qū)而言，各種重金屬對整個(gè)城區(qū)的危害程度不同。綜合羅綜合指數(shù)圖表以及重金元素在各功能區(qū)的空間分布，對數(shù)據(jù)進(jìn)行主成分分析，每種重金屬元

14、素的綜合污染指數(shù)與每個(gè)功能區(qū)具有線性關(guān)系，可求出每種元素在整個(gè)區(qū)的主成分貢獻(xiàn)率：as cd cr cuhg ni pb zn 0.02100.02580.09590.25510.39880.04530.03760.1206對各重金屬元素在整個(gè)城區(qū)的貢獻(xiàn)率進(jìn)行排序：hg、cu、zn、cr、ni、pb、cd、as。經(jīng)貢獻(xiàn)率排序可得，在整個(gè)城區(qū)內(nèi)，重金屬hg、cu兩種元素對環(huán)境的危害程度較大，而cd、as元素對環(huán)境的危害程度較小。其次，對五個(gè)功能區(qū)而言，各種重金屬對于每個(gè)功能區(qū)的危害程度不同。由模型一，可繪制出內(nèi)梅羅綜合指數(shù)圖表：根據(jù)圖表，可以直觀的觀察出環(huán)境污染嚴(yán)重由哪種重金屬元素造成。在c1（生

15、活區(qū)）內(nèi)，cr、zn含量較多，這些元素主要由生活中的化妝品、染發(fā)劑、燃煤、皮革、塑料門窗、蓄電池等生活用品產(chǎn)生。在c2（工業(yè)區(qū)）內(nèi)，as、hg含量較多，這些元素主要由工廠排放的廢渣、廢氣、廢水進(jìn)入環(huán)境，也告訴我們在工業(yè)生產(chǎn)中應(yīng)注意“三廢”的先處理后排放產(chǎn)生。在c3（山區(qū)）內(nèi)，ni含量較多，這是因?yàn)楹噺U水灌溉、植物生長吸收土壤中的鎳、動(dòng)植物的殘?bào)w腐爛、巖石風(fēng)化等原因造成山區(qū)重金屬的積累。在c4（交通區(qū)）內(nèi)重金屬元素含量比重都較大，這些元素的主要由機(jī)械中所含的重金屬、交通工具尾氣排放及輪胎磨損產(chǎn)生的大量含重金屬的有害氣體和粉塵等產(chǎn)生。在c5（公園綠地區(qū)）內(nèi)，hg、zn含量較多，這些元素主要由人類

16、在公園留下的廢棄物、大型娛樂機(jī)械產(chǎn)生。 4.3.問題（3）模型三：歸一化數(shù)值分析與污染源濃度目標(biāo)函數(shù)求解模型 污染源位置的確定與污染物的濃度有關(guān)，根據(jù)污染物的擴(kuò)散方式可以知道，污染物可能是由高濃度向低濃度擴(kuò)散，也可能是高地勢向低地勢擴(kuò)散。因此污染源可能在（正常范圍內(nèi)）污染物濃度較高的地方，也可能在地勢較高的地方。 按照同種濃度不同重金屬對環(huán)境影響程度標(biāo)準(zhǔn)，將所給的數(shù)據(jù)進(jìn)行歸一化處理，可以得到矩陣a（見附錄）。從得到的數(shù)據(jù)畫直方圖，可以看出圖中有幾處污染物的含量明顯偏高，污染源的位置可能就在其附近，污染物含量偏高的點(diǎn)可能為特殊點(diǎn)（可能是人為原因?qū)е?，如有人向此處堆積污染物），將此類點(diǎn)看為特殊點(diǎn)，

17、對它們進(jìn)行剔除。將污染源大概化為三個(gè)且分布在三個(gè)地理位置，然后按坐標(biāo)位置將其劃分為三個(gè)地理區(qū)域數(shù)據(jù)直方圖對上圖中的樣本點(diǎn)9,182,257對應(yīng)的地理坐標(biāo)，由圖2分析，對這三個(gè)點(diǎn)地理位置分別分析可知：樣本點(diǎn)9附近的生活區(qū)高度密集，受人類活動(dòng)影響因素較大，重金屬元素cr、zn濃度呈現(xiàn)較高的現(xiàn)象；樣本點(diǎn)182附近交通區(qū)與生活區(qū)交叉分布，重金屬元素hg、as濃度呈現(xiàn)較高的現(xiàn)象；樣本點(diǎn)257附近的工業(yè)區(qū)分布集中，重金屬元素hg、cu濃度呈現(xiàn)較高的現(xiàn)象；這些地區(qū)可能是由附近的交通要道、生活垃圾、工廠排放“三廢”導(dǎo)致的重金屬污染物富集，致使重金屬污染物濃度較高。同時(shí)綜合考慮到上述三處地方地勢較低，也可能是山

18、上的重金屬元素通過雨水沖刷、巖石風(fēng)化、大氣沉降等原因?qū)е碌摹Ｒ虼?，可以分兩情況考慮。假如是上述第一種情況導(dǎo)致的，我們可根據(jù)函數(shù)最優(yōu)規(guī)劃對此地區(qū)的污染源進(jìn)行較精確判斷，跟據(jù)此地區(qū)的地貌特征和污染物綜合指數(shù)的關(guān)系進(jìn)行數(shù)據(jù)多元擬合建立方程。設(shè)分別用數(shù)據(jù)進(jìn)行擬合得到c(t)=f(x),c(t)=f(y),c(t)=f(z)。則可得到污染物濃度合坐標(biāo)的對應(yīng)方程c（t）=p（q1,q2,q3）;然后，我們可根據(jù)無約束條件目標(biāo)函數(shù)求解方法，用matlab中的fminsearch函數(shù)讓目標(biāo)函數(shù)在出事猜測點(diǎn)8點(diǎn)所對應(yīng)的地理位置附近求污染物濃度最大點(diǎn)，即求min-c(t)（也即目標(biāo)函數(shù)在點(diǎn)8對應(yīng)坐標(biāo)的位置附近取最

19、大值），當(dāng)c（t）最大時(shí)求出的點(diǎn)坐標(biāo)就是理想的污染源所在點(diǎn)。圖2上面我們大致確定了現(xiàn)在污染比較嚴(yán)重的區(qū)域的位置，通過對地貌的分析發(fā)現(xiàn)這些位置大部分在山腳下，而在高層的重金屬會(huì)隨著雨水的沖刷流向低處，在相對低洼的地方富集，使得這些地方的重金屬污染濃度較高，因此估計(jì)這些地區(qū)的污染可能是有山上的某個(gè)地方堆積了大量的重金屬元素引起的，即可確定污染源的位置。利用歸一化和數(shù)據(jù)，分別與x,y,海拔高度建立擬合模型；經(jīng)分析得到的方程如下： y1=1.027e-012x3-2.71e-008x2 -0.0001086x+1.375 y2=-7.061e-013x3+3.706e-008x2-0.000641x+

20、4.305y3=-9.117e-007x3+0.0004927x2-0.07837x+4.185再利用微分方程的基本知識就可確定污染源大致在（15600，15286，118）。4.4.問題（4）：對模型三進(jìn)行評價(jià)：優(yōu)點(diǎn)：該模型綜合考慮了城區(qū)的地貌、重金屬污染物的傳播特征與其分布狀況，并對污染源的位置進(jìn)行分類思考，利用本方案來確定污染源的位置更具科學(xué)性。 缺點(diǎn)：對城區(qū)土壤地質(zhì)環(huán)境進(jìn)行調(diào)查時(shí)，采樣點(diǎn)的數(shù)據(jù)是在某個(gè)時(shí)刻測量出的重金屬元素的濃度，該模型缺乏對時(shí)間因素的考慮。不能夠很好的反映出城市地質(zhì)環(huán)境的演變過程。為了對模型三進(jìn)行優(yōu)化，還應(yīng)按一定梯度時(shí)間，在相同采樣點(diǎn)進(jìn)行測量各種重金屬元素的濃度大小和

21、對地表特征進(jìn)行逼近。設(shè)為擴(kuò)散域d的地表曲面方程，其定義域?yàn)?剖分為由凸圖形構(gòu)成的有限多子域之和, 對所有的、，有: ， 如果，令 那么，必然有 而且 其中為的邊界, 為域的近似中子域與有公共點(diǎn)( l) 的子域邊界，為指標(biāo)集。.作曲面的分片k 次插值,使之滿足:其中， 為曲面的次數(shù)的插值曲面，而 為插值基函數(shù)，并且有如果，令 那么，有：= 其中，為節(jié)點(diǎn)總數(shù)，域共p個(gè)子域；顯然，當(dāng)max0時(shí)，有：。擴(kuò)散域的曲面坐標(biāo)設(shè)曲面的參數(shù)形式為 如果 于是，有即雙方11對應(yīng)。如果再設(shè)曲面上任一點(diǎn)的位置向量為; 那么,曲面方程可以由, 給出，而且向量是線性無關(guān)的。如果令其中。如果再記 并且，設(shè) 那么，可以取，顯

22、然 最后，選取:與及和與及切于點(diǎn)的曲面和并且以 作為坐標(biāo)曲面，而且有曲面坐標(biāo)系擴(kuò)散方程：設(shè)為擴(kuò)散域內(nèi)的，曲面坐標(biāo)擴(kuò)散系方程和濃度的關(guān)系，同樣通過matlab中的fminsearch函數(shù)，可不同時(shí)刻求出理想的污染源的位置。由污染源位置的變化和各種重金屬元素對環(huán)境影響程度的變化可以看出人類活動(dòng)對環(huán)境影響程度的變化。五、模型的優(yōu)缺點(diǎn)優(yōu)點(diǎn)：（1）、模型中不僅考慮了不同元素對環(huán)境的影響程度，并且考慮到不同功能區(qū)的環(huán)境污染程度，將不同功能的污染程度也進(jìn)行了排序。（2）、對于重金屬污染的主要原因，先整體進(jìn)行分析，后劃分五個(gè)模塊對五個(gè)功能區(qū)不同元素對其不同影響程度進(jìn)行了分析。（3）、綜合考慮了城區(qū)的地貌、重金

23、屬污染物的傳播特征與其分布狀況，并對污染源的位置進(jìn)行分類思考，利用本方案來確定污染源的位置更具科學(xué)性。缺點(diǎn)：（1）、重金屬元素空間分布圖中，缺乏海拔高度變量對其分布的影響。（2）、對污染源空間位置確定時(shí)，求得的坐標(biāo)為該區(qū)域內(nèi)較理想的空間坐標(biāo)，而非污染源的準(zhǔn)確位置。 參考文獻(xiàn)1 韓中庚，數(shù)學(xué)建模方法及其應(yīng)用m，北京：高等教育出版社，2005年。2 姜啟源，謝金星，葉俊，數(shù)學(xué)模型，北京：高等教育出版社，2004年。3 宋葉志，賈東永，matlab數(shù)值分析與應(yīng)用，北京：機(jī)械工業(yè)出版社，20094 劉衍君，湯慶新，白振華，張秀玲，張保華，基于地質(zhì)累積與內(nèi)梅羅指數(shù)的耕地重金屬污染研究，附錄：程序代碼1：

24、x=shuju(:,2); y=shuju(:,3); z=shuju(:,6); for i=1:319if shuju(i,5)=1plot3(x(i),y(i),z(i),*)else if shuju(i,5)=2plot3(x(i),y(i),z(i),+)else if shuju(i,5)=3plot3(x(i),y(i),z(i),x)else if shuju(i,5)=4plot3(x(i),y(i),z(i),o) else plot3(x(i),y(i),z(i),p)endhold onendendendendhold onx0=shuju(:,2); y0=shuj

25、u(:,3); z0=shuju(:,6); x=0:100:30000; y=0:100:20000;x,y=meshgrid(x,y); z=griddata(x0,y0,z0,x,y,cubic);subplot(1,2,1)mesh(x,y,z),subplot(1,2,2)cc,hh=contour(x,y,z); clabel(cc,hh)% c,h = contour(peaks); clabel(c,h), colorbarhold onx=shuju(:,2); y=shuju(:,3); z=shuju(:,6); for i=1:319if shuju(i,5)=1plo

26、t(x(i),y(i),*)else if shuju(i,5)=2plot(x(i),y(i),+)else if shuju(i,5)=3plot(x(i),y(i),x)else if shuju(i,5)=4plot(x(i),y(i),o)else plot(x(i),y(i),p)endhold onendendendend程序代碼2：b1=3.6;b2=130;b3=31;b4=13.2;b5=35;b6=12.3;b7=31;b8=69;i11=log2(a1(:,1)./(1.5*b1);i21=log2(a1(:,2)./(1.5*b2);i31=log2(a1(:,3).

27、/(1.5*b3);i41=log2(a1(:,4)./(1.5*b4);i51=log2(a1(:,5)./(1.5*b5);i61=log2(a1(:,6)./(1.5*b6);i71=log2(a1(:,7)./(1.5*b7);i81=log2(a1(:,8)./(1.5*b8);b1=i11 i21 i31 i41 i51 i61 i71 i81;c11=max(i11);d11=mean(i11);c21=max(i21);d21=mean(i21);c31=max(i31);d31=mean(i31);c41=max(i41);d41=mean(i41);c51=max(i51)

28、;d51=mean(i51);c61=max(i61);d61=mean(i61);c71=max(i71);d71=mean(i71);c81=max(i81);d81=mean(i81);p11=(c112+d112)/20.5;p21=(c212+d212)/20.5;p31=(c312+d312)/20.5;p41=(c412+d412)/20.5;p51=(c512+d512)/20.5;p61=(c612+d612)/20.5;p71=(c712+d712)/20.5;p81=(c812+d812)/20.5;c1=p11 p21 p31 p41 p51 p61 p71 p81i1

29、2=log2(a2(:,1)./(1.5*b1);i22=log2(a2(:,2)/(1.5*b2);i32=log2(a2(:,3)./(1.5*b3);i42=log2(a2(:,4)./(1.5*b4);i52=log2(a2(:,5)./(1.5*b5);i62=log2(a2(:,6)./(1.5*b6);i72=log2(a2(:,7)./(1.5*b7);i82=log2(a2(:,8)./(1.5*b8);b2=i12 i22 i32 i42 i52 i62 i72 i82;c12=max(i12);d12=mean(i12);c22=max(i22);d22=mean(i22

30、);c32=max(i32);d32=mean(i32);c42=max(i42);d42=mean(i42);c52=max(i52);d52=mean(i52);c62=max(i62);d62=mean(i62);c72=max(i72);d72=mean(i72);c82=max(i82);d82=mean(i82);p12=(c122+d122)/20.5;p22=(c222+d222)/20.5;p32=(c322+d322)/20.5;p42=(c422+d422)/20.5;p52=(c522+d522)/20.5;p62=(c622+d622)/20.5;p72=(c722+

31、d722)/20.5;p82=(c822+d822)/20.5;c2=p12 p22 p32 p42 p52 p62 p72 p82i13=log2(a3(:,1)./(1.5*b1);i23=log2(a3(:,2)/(1.5*b2);i33=log2(a3(:,3)./(1.5*b3);i43=log2(a3(:,4)./(1.5*b4);i53=log2(a3(:,5)./(1.5*b5);i63=log2(a3(:,6)./(1.5*b6);i73=log2(a3(:,7)./(1.5*b7);i83=log2(a3(:,8)./(1.5*b8);b3=i13 i23 i33 i43

32、i53 i63 i73 i83;c13=max(i13);d13=mean(i13);c23=max(i23);d23=mean(i23);c33=max(i33);d33=mean(i33);c43=max(i43);d43=mean(i43);c53=max(i53);d53=mean(i53);c63=max(i63);d63=mean(i63);c73=max(i73);d73=mean(i73);c83=max(i83);d83=mean(i83);p13=(c132+d132)/20.5;p23=(c232+d232)/20.5;p33=(c332+d332)/20.5;p43=(

33、c432+d432)/20.5;p53=(c532+d532)/20.5;p63=(c632+d632)/20.5;p73=(c732+d732)/20.5;p83=(c832+d832)/20.5;c3=p13 p23 p33 p43 p53 p63 p73 p83i14=log2(a4(:,1)./(1.5*b1);i24=log2(a4(:,2)/(1.5*b2);i34=log2(a4(:,3)./(1.5*b3);i44=log2(a4(:,4)./(1.5*b4);i54=log2(a4(:,5)./(1.5*b5);i64=log2(a4(:,6)./(1.5*b6);i74=l

34、og2(a4(:,7)./(1.5*b7);i84=log2(a4(:,8)./(1.5*b8);b4=i14 i24 i34 i44 i54 i64 i74 i84;c14=max(i14);d14=mean(i14);c24=max(i24);d24=mean(i24);c34=max(i34);d34=mean(i34);c44=max(i44);d44=mean(i44);c54=max(i54);d54=mean(i54);c64=max(i64);d64=mean(i64);c74=max(i74);d74=mean(i74);c84=max(i84);d84=mean(i84);

35、p14=(c142+d142)/20.5;p24=(c242+d242)/20.5;p34=(c342+d342)/20.5;p44=(c442+d442)/20.5;p54=(c542+d542)/20.5;p64=(c642+d642)/20.5;p74=(c742+d742)/20.5;p84=(c842+d842)/20.5;c4=p14 p24 p34 p44 p54 p64 p74 p84i15=log2(a5(:,1)./(1.5*b1);i25=log2(a5(:,2)/(1.5*b2);i35=log2(a5(:,3)./(1.5*b3);i45=log2(a5(:,4)./

36、(1.5*b4);i55=log2(a5(:,5)./(1.5*b5);i65=log2(a5(:,6)./(1.5*b6);i75=log2(a5(:,7)./(1.5*b7);i85=log2(a5(:,8)./(1.5*b8);b5=i15 i25 i35 i45 i55 i65 i75 i85;c15=max(i15);d15=mean(i15);c25=max(i25);d25=mean(i25);c35=max(i35);d35=mean(i35);c45=max(i45);d45=mean(i45);c55=max(i55);d55=mean(i55);c65=max(i65);

37、d65=mean(i65);c75=max(i75);d75=mean(i75);c85=max(i85);d85=mean(i85);p15=(c152+d152)/20.5;p25=(c252+d252)/20.5;p35=(c352+d352)/20.5;p45=(c452+d452)/20.5;p55=(c552+d552)/20.5;p65=(c652+d652)/20.5;p75=(c752+d752)/20.5;p85=(c82+d82)/20.5;c5=p15 p25 p35 p45 p55 p65 p75 p85a1=9.17 287.00 43.94 45.77 168.0

38、0 19.70 62.74 223.16 5.72 193.70 80.35 26.57 111.00 19.80 57.64 89.08 11.45 1044.50 94.78 136.97 202.00 22.30 472.48 602.04 7.84 347.90 57.65 97.14 213.00 19.60 70.82 307.24 8.50 614.00 744.46 130.55 156.00 32.80 228.64 1013.47 5.51 257.20 54.64 29.01 104.00 13.20 87.68 223.27 9.39 325.80 172.29 104

39、.89 82.00 31.50 90.90 429.29 4.09 90.50 35.02 11.82 16.00 10.40 29.09 46.84 6.35 532.00 57.51 83.76 191.00 19.50 73.46 297.14 3.50 396.30 138.37 58.97 170.00 24.20 91.76 2893.47 4.29 449.10 67.22 51.64 315.00 15.70 106.97 294.69 7.41 337.30 77.27 248.85 90.00 20.10 99.58 210.00 2.91 265.00 35.66 29.

40、39 24.64 9.23 60.54 122.96 3.30 488.00 51.18 34.55 37.50 10.80 54.62 125.92 6.14 227.00 42.15 67.04 49.29 16.31 34.28 82.96 8.06 113.10 52.40 20.81 65.36 19.69 29.56 62.24 4.69 132.80 36.11 17.28 20.36 15.53 37.03 82.76 2.34 87.60 18.46 9.73 13.93 8.89 24.43 43.37 6.56 245.60 36.73 61.30 55.71 14.18

41、 47.24 114.29 8.23 756.40 42.73 87.52 63.00 19.26 88.74 184.69 8.90 307.30 54.39 57.21 326.00 25.72 131.93 256.94 5.41 178.90 29.54 23.73 52.00 9.89 49.84 118.88 7.78 315.50 49.76 28.03 550.00 18.95 45.73 109.29 5.62 134.60 25.33 19.10 45.00 11.66 40.50 87.14 6.05 365.00 35.92 30.91 110.00 17.91 44.

42、75 147.35 4.17 310.70 40.70 42.64 58.00 14.99 115.59 177.76 6.26 387.20 38.03 30.06 85.00 21.97 54.98 142.55 5.00 196.50 50.03 18.56 29.00 25.82 33.90 84.80 4.58 129.00 31.09 18.93 38.00 15.10 29.76 69.80 5.41 204.90 40.16 25.86 41.00 16.24 33.28 80.00 7.56 206.60 55.79 24.44 22.00 28.63 29.73 81.73

43、 5.41 250.10 34.32 16.25 15.00 18.53 41.29 90.20 5.83 89.20 54.90 23.28 19.00 21.76 24.95 53.57 5.20 334.30 47.05 40.37 76.00 18.40 61.03 211.47 5.20 113.70 41.85 20.73 12.00 19.00 31.87 57.27 4.38 246.30 29.90 21.84 84.00 14.10 59.00 95.17 5.41 158.00 46.86 24.02 31.00 19.30 36.27 79.24 8.67 245.70

44、 47.79 27.55 35.00 18.40 53.42 98.81 6.47 86.80 41.12 15.46 23.00 15.90 37.53 70.18 7.12 367.80 92.02 49.80 97.00 16.30 41.26 321.12 10.97 248.50 40.61 61.52 81.00 17.20 76.19 168.05 9.81 171.80 75.38 163.20 30.00 26.30 45.27 125.16 2.77 236.20 42.67 16.35 62.00 9.40 41.88 149.52 6.47 282.90 52.68 2

45、0.34 25.00 22.90 32.53 103.50 a2=6.56 223.90 40.08 25.17 950.00 15.40 32.28 117.35 14.08 1092.90 67.96 308.61 1040.00 28.20 434.80 966.73 9.62 1066.20 285.58 2528.48 13500.00 41.70 381.64 1417.86 21.87 424.50 73.40 59.72 1520.00 27.80 83.70 175.71 18.38 630.00 96.68 114.81 645.00 34.80 130.36 1626.0

46、2 10.53 635.30 64.03 101.35 190.00 28.30 162.64 615.10 2.34 353.00 24.53 12.70 11.79 9.00 58.80 89.08 9.35 407.50 55.54 61.83 112.00 24.05 66.82 208.27 4.79 289.70 42.12 41.10 57.00 17.49 76.80 237.45 1.61 295.80 15.40 18.35 19.00 4.27 40.42 106.53 4.79 178.60 32.31 29.01 45.00 14.26 47.98 104.90 8.

47、67 423.30 61.41 117.83 216.00 30.30 178.88 293.16 7.12 967.70 36.41 70.71 66.00 16.87 119.35 457.96 4.58 445.70 28.32 132.05 139.00 12.91 87.44 303.06 8.23 420.60 35.81 73.86 561.00 22.49 60.91 241.63 4.58 257.40 42.25 44.84 260.00 19.58 46.90 192.65 8.23 350.30 59.07 36.58 177.00 24.57 52.56 170.00

48、 3.97 306.20 38.30 36.22 115.00 13.54 60.63 145.51 10.74 306.50 53.98 97.05 255.00 27.91 69.40 177.86 3.77 515.60 28.74 77.47 1801.00 15.10 70.18 178.57 10.27 175.90 57.88 37.50 65.00 26.97 45.34 100.41 8.90 117.50 32.54 20.51 41.00 14.47 41.45 56.33 3.37 379.10 32.13 37.81 49.00 13.74 70.37 208.37

49、6.69 320.30 87.90 66.99 78.00 21.90 123.60 230.11 5.00 551.10 71.75 71.23 74.00 22.30 111.62 209.50 9.58 201.00 39.11 25.26 125.00 10.70 99.35 90.07 3.17 203.80 26.57 24.29 23.00 11.90 41.69 94.34 3.77 355.20 26.25 32.87 14.00 9.60 59.22 117.87 3.37 240.50 28.04 22.63 261.00 11.70 35.74 92.46 6.05 1

50、93.00 40.60 24.88 27.00 14.40 33.53 84.86 3.77 467.10 49.03 34.44 45.00 15.40 60.83 132.86 5.41 364.70 40.34 40.93 79.00 18.70 83.32 175.34 8.23 409.90 44.67 66.92 80.00 36.00 96.85 197.63 5.41 302.50 34.22 27.60 408.00 14.80 68.70 218.24 7.78 114.50 56.38 26.96 36.00 22.40 31.24 75.91 6.47 165.20 73.40 42.73 40.00 19.70 84.13 95.694.09 127.00 27.58 23.99 30.00 11.93 57.47 85.61 2.72 70.90 19.45 9.12 15.00 7.09 22.73 32.86 3.69 150.70 59.61 19.00 34.29 24.98 38.47 89.08 6.14 100.30 37.49 20.23 34.29 14.85 29.29 61.94 10.99 109.80 56.07 69.06 58.93 20.70 38.87 63.27 7.41 114.30 4