承 諾 書

我們仔細(xì)閱讀了中國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽的競(jìng)賽規(guī)則.

我們完全明白，在競(jìng)賽開始后參賽隊(duì)員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊(duì)外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問(wèn)題。

我們知道，抄襲別人的成果是違反競(jìng)賽規(guī)則的, 如果引用別人的成果或其他公開的資料（包括網(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)填寫）： 我們的參賽報(bào)名號(hào)為（如果賽區(qū)設(shè)置報(bào)名號(hào)的話）： 所屬學(xué)校（請(qǐng)?zhí)顚懲暾娜?

參賽隊(duì)員 (打印并簽名) ：1.

2.

3. 指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人 (打印并簽名)：

日期： 年 月 日

賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：

1

2011高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽

編 號(hào) 專 用 頁(yè)

評(píng) 閱 人 評(píng) 分 備 注 賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：

賽區(qū)評(píng)閱記錄（可供賽區(qū)評(píng)閱時(shí)使用）：

全國(guó)統(tǒng)一編號(hào)（由賽區(qū)組委會(huì)送交全國(guó)前編號(hào)）：

全國(guó)評(píng)閱編號(hào)（由全國(guó)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：

2

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

摘要

本題主要通過(guò)對(duì)重金屬污染濃度的分析，畫出濃度的空間分布圖，建立相關(guān)的模型，得出該地區(qū)的重金屬污染程度、污染原因以及污染源的所在位置。

針對(duì)不同重金屬具有不同的污染源，對(duì)不同區(qū)域重金屬元素濃度值進(jìn)行分析。運(yùn)用克拉格差值法，畫出各重金屬的濃度等高線圖。從水平和垂直方向考慮濃度分布，在水平面與海拔上對(duì)金屬濃度分別分析，得出在三維空間內(nèi)濃度的分布狀況?？紤]到區(qū)域內(nèi)各種重金屬元素污染程度，及不同區(qū)域內(nèi)元素的綜合污染狀況不同，先運(yùn)用地積累法，對(duì)各區(qū)域內(nèi)不同元素分別考慮其污染程度，再對(duì)每個(gè)區(qū)域用內(nèi)梅羅綜合污染指數(shù)法得出該地區(qū)綜合污染程度。

根據(jù)給出的樣本數(shù)據(jù)對(duì)各種重金屬元素之間的相關(guān)性進(jìn)行分析，得出相關(guān)性矩陣，對(duì)于相關(guān)性比較明顯的重金屬，我們可以認(rèn)為其來(lái)源大致相同，即有相同的污染源。然后采用KMO法和Bartlett法對(duì)原始數(shù)據(jù)進(jìn)行主成分分析適宜性檢驗(yàn)，在此基礎(chǔ)上進(jìn)一步對(duì)樣本數(shù)據(jù)進(jìn)行主成分分析，確定影響比較顯著的幾種重金屬元素，從而簡(jiǎn)化模型，方便分析重金屬污染的主要原因。

在分析重金屬污染傳播的特征時(shí)，可以結(jié)合大氣擴(kuò)散模型和水體對(duì)流擴(kuò)散模型的偏微分方程建立重金屬元素的對(duì)流擴(kuò)散模型。為求具體污染源的位置，在對(duì)擴(kuò)散模型簡(jiǎn)化的基礎(chǔ)上，利用反演擴(kuò)散模型方程來(lái)求污染源。反演方程可以用遺傳算法來(lái)實(shí)現(xiàn)，把問(wèn)題轉(zhuǎn)化為建立一個(gè)非線性最優(yōu)化問(wèn)題。

比較用差分方法反演擴(kuò)散模型與遺傳算法反演擴(kuò)散模型的優(yōu)缺點(diǎn)，得出差分方法是建立在把偏微分方程轉(zhuǎn)化為差分方程進(jìn)而求解各個(gè)點(diǎn)處的濃度隨時(shí)間變化差異性確定污染源，針對(duì)樣本數(shù)據(jù)分布不規(guī)律，需要進(jìn)行插值處理，會(huì)給模型的精度帶來(lái)較大影響。而遺傳算法具有極好的全局優(yōu)化性能,減少了陷入局部極值的風(fēng)險(xiǎn), 特別對(duì)非線性優(yōu)化問(wèn)題具有很強(qiáng)的適應(yīng)性，所以采用遺傳算法反演擴(kuò)散方程，具有較好的合理性。但由于未考慮模型里的對(duì)流項(xiàng)以及海拔的影響，導(dǎo)致其誤差較大，如能提供該地區(qū)具體的工業(yè)排放區(qū)及交通要道等的具體密集情況，進(jìn)一步減小污染源位置的范圍和個(gè)數(shù)的初始值，可以提高遺傳算法得出的污染源的具體位置的準(zhǔn)確性。

關(guān)鍵詞：等高線圖 相關(guān)性 主成分分析 遺傳算法 反演 擴(kuò)散方程

3

目錄

1、 問(wèn)題重述……………………………………………………………….3 2、 問(wèn)題分析……………………………………………………………….3

2.1、對(duì)問(wèn)題一的分析…………………………………………………3 2.2、對(duì)問(wèn)題二的分析…………………………………………………3 2.3、對(duì)問(wèn)題三的分析…………………………………………………4 2.4、對(duì)問(wèn)題四的分析…………………………………………………4 3、 模型假設(shè)……………………………………………………………….4 4、 符號(hào)說(shuō)明……………………………………………………………….4 5、 模型的建立與求解…………………………………………………….5

5.1、問(wèn)題一的模型建立與求解………………………………………5 5.1.1、重金屬空間分布…………………………………………...5 5.1.2、重金屬的污染程度的分析………………………………...8 5.2、問(wèn)題二的模型建立和求解………………………………………11 5.2.1、重金屬相關(guān)性分析………………………………………...11 5.2.2、土壤重金屬主成分分析……………………………………13 5.3、問(wèn)題三的模型建立和求解………………………………………15

5.3.1 污染物傳播特征的分析……………………………………15

5.3.2模型建立和求解……………………………………………15 5.3.3模型可行性分析……………………………………………18 5.4問(wèn)題四的分析和求解………………………………………………19

5.4.1優(yōu)缺點(diǎn)分析…………………………………………………19 5.4.2模型優(yōu)化……………………………………………………19

6、參考文獻(xiàn)……………………………………………………………… 20 7、附錄…………………………………………………………………….21

4

1、問(wèn)題重述

隨著城市經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷增加，人類活動(dòng)對(duì)城市環(huán)境質(zhì)量的影響日顯突出。對(duì)城市土壤地質(zhì)環(huán)境異常的查證，以及如何應(yīng)用查證獲得的海量數(shù)據(jù)資料開展城市環(huán)境質(zhì)量評(píng)價(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)影響的程度不同。

現(xiàn)對(duì)某城市城區(qū)土壤地質(zhì)環(huán)境進(jìn)行調(diào)查。為此，將所考察的城區(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列出了采樣點(diǎn)的位置、海拔高度及其所屬功能區(qū)等信息，附件2列出了8種主要重金屬元素在采樣點(diǎn)處的濃度，附件3列出了8種主要重金屬元素的背景值。

現(xiàn)要求你們通過(guò)數(shù)學(xué)建模來(lái)完成以下任務(wù)：

(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)題？

2、問(wèn)題分析

2.1 對(duì)問(wèn)題一的分析

題目要求統(tǒng)計(jì)出8重金屬元素在城區(qū)的空間分布，需要考慮重金屬元素在X，Y，Z三維空間的濃度分布特征。由附件給出的數(shù)據(jù)分析得，當(dāng)X，Y的坐標(biāo)確定之后，該點(diǎn)的海拔高度Z相應(yīng)確定下來(lái)。所以在統(tǒng)計(jì)重金屬的空間分布時(shí)，可以先考慮元素濃度和平面位置的關(guān)系，再通過(guò)海拔分布圖分析元素分布與海拔之間的關(guān)系。

對(duì)這八種重金屬元素分別單獨(dú)考慮，利用附件中各元素在城區(qū)的樣本濃度，采用克立格插值法，應(yīng)用SUPFER 8.0軟件對(duì)各元素含量分布進(jìn)行空間分析，可分別繪制出8種重金屬元素的空間等高線分布圖。

利用附件給出了的各重金屬元素在當(dāng)?shù)氐谋尘爸?，可以運(yùn)用地積累指數(shù)法來(lái)分析該城區(qū)不同區(qū)域8種重金屬元素各自的污染程度。同時(shí)，也可以運(yùn)用內(nèi)梅羅指數(shù)法綜合評(píng)判該城區(qū)各區(qū)域內(nèi)的重金屬污染程度。

2.2 對(duì)問(wèn)題二的分析

根據(jù)附件給出的數(shù)據(jù)，我們分區(qū)域考慮，確定每個(gè)區(qū)域里8種重金屬元素相互之間的相關(guān)性，相關(guān)性較顯著的重金屬元素可以認(rèn)為來(lái)自于同一個(gè)污染源，由

5

此可以大致確定每個(gè)區(qū)域污染源的種類。因?yàn)槊總€(gè)區(qū)域中重金屬的種類較多，相互之間又具有比較明顯的相關(guān)性，所以我們可以在采用KMO法和Bartlett法對(duì)原始數(shù)據(jù)進(jìn)行主成分分析適宜性檢驗(yàn)的基礎(chǔ)上，進(jìn)一步考慮使用主成分分析法，由此分析出每個(gè)區(qū)域影響比較顯著的若干種重金屬，從而確定出重金屬污染的主要原因。

2.3 對(duì)問(wèn)題三的分析

重金屬元素傳播途徑主要有水流傳播、土壤滲透、大氣擴(kuò)散等。在城區(qū)我們主要考慮以大氣擴(kuò)散方式傳播。由此，我們通過(guò)對(duì)流擴(kuò)散方程建立模型，從中體現(xiàn)重金屬傳播的主要特性。在隨時(shí)間的傳播過(guò)程中，污染源所在地的濃度隨時(shí)間的變化率較大。由此可以通過(guò)對(duì)模型的反演確定污染源的大概位置。

2.4 對(duì)問(wèn)題四的分析

前三個(gè)問(wèn)題的求解模型雖然可以比較合理地解決問(wèn)題，但是也存在一定的局限性。由此我們可以進(jìn)一步收集該城區(qū)土壤成分的信息，建立更優(yōu)化的模型。

3、模型假設(shè)

1、固定時(shí)間內(nèi)污染元素濃度保持穩(wěn)定。

2、元素污染濃度在水平面與海拔高度的分布成的相關(guān)性。 3、污染元素的擴(kuò)散只受水平面的擴(kuò)散影響。

4、污染元素在水平面的擴(kuò)散在各方向的擴(kuò)散速度相同。 5、污染源的位置大致固定，不隨時(shí)間變化。

4、符號(hào)說(shuō)明

Ci： 重金屬污染物i的實(shí)測(cè)濃度 Bi： 重金屬元素n的背景值 Igeo：地積累指數(shù)

PI： 重金屬污染物i的污染指數(shù) Si： 重金屬污染物i的國(guó)家標(biāo)準(zhǔn)值

Pmax：某區(qū)域中重金屬單項(xiàng)污染指數(shù)最大值 Pave：某區(qū)域所有重金屬單項(xiàng)污染指數(shù)的平均值 IPI： 某區(qū)域土壤中重金屬綜合污染指數(shù)

?x： x方向上的網(wǎng)格步長(zhǎng)

Ex ： x方向的擴(kuò)散系數(shù) Ey ： x方向的擴(kuò)散系數(shù) Ez ： x方向的擴(kuò)散系數(shù)

6

5、模型的建立和求解

5.1問(wèn)題一的模型建立和求解 5.1.1 重金屬空間分布的求解

克里格插值法（Kriging）是用協(xié)方差函數(shù)和變異函數(shù)來(lái)確定高程變量隨空間距離而變化的規(guī)律,以距離為自變量的變異函數(shù),計(jì)算相鄰高程值關(guān)系權(quán)值,在有限區(qū)域內(nèi)對(duì)區(qū)域化變量進(jìn)行無(wú)偏最優(yōu)估計(jì)的一種方法，是地統(tǒng)計(jì)學(xué)的主要方法之一。

根據(jù)采樣點(diǎn)各重金屬元素的含量分布，采用克立格插值法，應(yīng)用SUPFER 8.0軟件對(duì)海拔高度和8種重金屬的含量分布進(jìn)行空間分析，可繪出以下空間分布圖。

圖5-1-1 海拔分布圖

圖5-1-2 As空間分布圖

由圖5-1-2分析得，As元素濃度較高部分在左圖呈明顯黑色的三個(gè)區(qū)域中，說(shuō)明這三個(gè)區(qū)域As元素污染較嚴(yán)重。同時(shí)考慮濃度與海拔的關(guān)系，右圖呈現(xiàn)各濃度主要分布在海拔高度為0—100m之間的范圍內(nèi)，且與海拔無(wú)特別關(guān)系。

7

圖5-1-3 Cd空間分布圖

由圖5-1-3分析得，Cd元素濃度較高部分在左圖中的幾個(gè)顏色較深的區(qū)域，且該元素在城區(qū)內(nèi)的污染都比較嚴(yán)重。由右圖得在海拔高度為0—100m的范圍內(nèi)均有濃度分布，且濃度較高部分分布在海拔較低區(qū)域內(nèi)。

圖5-1-4 Cr空間分布圖

由圖5-1-4分析得，Cr元素在城區(qū)內(nèi)總體分布較少，但有一塊很明顯的黑色區(qū)域，表明污染主要集中在這一塊。同時(shí)，右圖的海拔圖顯示污染主要分布于0-100m海拔之間，且濃度量大致在0-100 (μg/g)之間。

圖5-1-5 Cu空間分布圖

由圖5-1-5分析得，Cu元素在城區(qū)內(nèi)的總體分布較少，高濃度主要集中于兩塊明顯黑色的區(qū)域，表明這些區(qū)域污染較嚴(yán)重。右圖的海拔圖顯示污染主要分布于0-100m海拔之間，且與海拔無(wú)特別關(guān)系。

8

圖5-1-6 Hg空間分布

由圖5-1-6分析得，Hg元素在城區(qū)內(nèi)的總體分布較少，高濃度主要集中于三塊明顯黑色的區(qū)域，表明這些區(qū)域污染較嚴(yán)重。右圖的海拔圖顯示污染主要分布于0-100m海拔之間，且與海拔無(wú)特別關(guān)系。

圖5-1-7 Ni空間分布圖

由圖5-1-7分析得，Ni元素在城區(qū)內(nèi)的總體分布較少，高濃度主要集中于二塊明顯黑色的區(qū)域，表明這些區(qū)域污染較嚴(yán)重。右圖的海拔圖顯示污染主要分布于0-100m海拔之間， 且與海拔無(wú)特殊關(guān)系。

圖5-1-8 Pb空間分布圖

由圖5-1-8分析得，Pb元素濃度較高部分在左圖中的兩個(gè)顏色較深的區(qū)域，說(shuō)明這兩個(gè)區(qū)域As元素污染較嚴(yán)重。由右圖得在海拔高度為0—100m的范圍內(nèi)均有濃度分布，且濃度較高部分大致分布在海拔較低區(qū)域內(nèi)。

9

圖5-1-9 Zn空間分布圖

由圖5-1-9分析得，Zn元素濃度較高部分在左圖中的幾個(gè)顏色較深的區(qū)域，說(shuō)明這幾個(gè)區(qū)域Zn元素污染較嚴(yán)重。由右圖得在海拔高度為0—100m的范圍內(nèi)均有濃度分布，且與海拔無(wú)特殊關(guān)系。

綜上所述，因?yàn)楦鞣N重金屬有不同的來(lái)源，比如交通污染，工業(yè)污染，生活垃圾污染等，所以呈現(xiàn)出不同的區(qū)域分布。但從海拔方面來(lái)分析，其主要分布在0—100m的低海拔區(qū)域。

5.1.2、重金屬的污染程度的分析 1）模型建立與求解

我們使用地積累指數(shù)法建立模型。Muller曾提出地積累指數(shù)定量評(píng)價(jià)沉積物中的重金屬污染程度，并規(guī)定了相應(yīng)的污染程度級(jí)別劃分標(biāo)準(zhǔn)（表5-1-1），這種評(píng)價(jià)方法也可以用來(lái)評(píng)價(jià)土壤中的重金屬的污染程度及其分級(jí)情況。其計(jì)算公式為：

Igeo?log2?C/KB?

ii式中Ci為重金屬元素i在樣本中的濃度，Bi為重金屬元素i的背景值，K是為考慮各地巖石差異可能會(huì)引起背景值的變動(dòng)而取的系數(shù)（一般取值K=1.5）。

表5-1-1 Muller 地積累指數(shù)分級(jí)表

根據(jù)附件表格給出的各重金屬元素在多個(gè)采樣點(diǎn)的濃度，同時(shí)我們也知道各個(gè)采樣點(diǎn)對(duì)應(yīng)的區(qū)域，所以我們可以用Matlab計(jì)算出每種重金屬在每個(gè)區(qū)域的平均濃度Ci。Bi取各元素背景值的平均值代入計(jì)算。則可以算出各元素在各個(gè)區(qū)域的地積累指數(shù)Igeo。

10

表5-1-2生活區(qū)重金屬污染程度

表5-1-3工業(yè)區(qū)重金屬污染程度

表5-1-4山區(qū)重金屬污染程度

表5-1-5主干道路區(qū)重金屬污染程度

11

表5-1-6公園綠地區(qū)重金屬污染程度

2）為了全面反映各污染物對(duì)土壤的不同作用，突出高濃度污染物對(duì)環(huán)境質(zhì)量的影響，采用內(nèi)梅羅綜合污染指數(shù)法再次建立一個(gè)模型。內(nèi)梅羅指數(shù)不僅考慮到各種污染物的平均污染狀況，并且強(qiáng)調(diào)最嚴(yán)重的污染物對(duì)污染指數(shù)的影響。因此，可采用單因子污染指數(shù)( PI) 和內(nèi)梅羅綜合污染指數(shù)( IPI) 對(duì)城區(qū)土壤重金屬的污染程度進(jìn)行綜合評(píng)價(jià)。

單因子污染指數(shù)法：

PI=Ci/Si

式中，PI 為重金屬污染物i的污染指數(shù);Ci為重金屬污染物i的實(shí)測(cè)濃度，Si為重金屬污染物i的國(guó)家標(biāo)準(zhǔn)值。

內(nèi)梅羅綜合污染指數(shù)法:

22Pmax?PaveIPI= 2式中，Pmax為一個(gè)區(qū)域采樣點(diǎn)中重金屬單項(xiàng)污染指數(shù)最大值；Pave為一個(gè)區(qū)域所有重金屬單項(xiàng)污染指數(shù)的平均值。IPI為一個(gè)區(qū)域土壤中重金屬綜合污染指數(shù)。評(píng)判標(biāo)準(zhǔn)如下：IPI≤1，表示無(wú)污染；1 < IPI≤2，表示輕度污染；2 < IPI≤3，表示中度污染；IPI > 3，表示嚴(yán)重污染。

由此可以進(jìn)一步計(jì)算出內(nèi)梅羅綜合污染指數(shù)( IPI)。我們把城區(qū)分為5個(gè)區(qū)域分別計(jì)算。計(jì)算過(guò)程中Ci取一個(gè)區(qū)域中重金屬i的采樣值的平均濃度，Si取為重金屬i的國(guó)家標(biāo)準(zhǔn)值。

通過(guò)Matlab編程計(jì)算出的結(jié)果如下表5-1-7所示。

表5-1-7各區(qū)域綜合污染程度

12

5.2、 問(wèn)題二的模型建立和求解 5.2.1、重金屬元素相關(guān)性分析

城市土壤重金屬來(lái)源于成土母質(zhì)和人類活動(dòng)，同一來(lái)源的重金屬之間存在著相關(guān)性，根據(jù)相關(guān)性可以判斷土壤重金屬污染來(lái)源是否相同。如果重金屬之間存在顯著的正相關(guān)，則其來(lái)源可能相同，否則來(lái)源可能不止一個(gè)。

應(yīng)用Matlab編程，可以計(jì)算出各個(gè)區(qū)域內(nèi)8種重金屬相互之間的相關(guān)性，并對(duì)其相關(guān)性做出分析。

表5-2-1生活區(qū)重金屬元素相關(guān)性

在生活區(qū)采集的土壤樣品中，As和Ni，Cd和Pb呈現(xiàn)比較明顯的相關(guān)性，相關(guān)系數(shù)分別是0.605和0.802.，表明As和Ni可能具有相同的污染源，Cd和Pb可能具有相同的污染源。但其他元素之間的相關(guān)性都不明顯，所以生活區(qū)的重金屬來(lái)源相對(duì)復(fù)雜，很難于某種污染源進(jìn)行確認(rèn)。

表5-2-2工業(yè)區(qū)重金屬元素相關(guān)性

在工業(yè)區(qū)采集的土壤樣品中，Cd、Cr、Cu、Hg、Pb、Zn兩兩之間都有比較明顯的相關(guān)性，表明這六種重金屬元素可能具有相同的污染源，As和Ni具有明顯的相關(guān)性，可能具有相同的污染源。

13

表5-2-3山區(qū)重金屬元素相關(guān)性

在山區(qū)采集的土壤樣品中，Cd、Pb、Zn兩兩之間呈現(xiàn)明顯的相關(guān)性，所以此三者可能有相同的污染源。As和Cu、Cu和Hg之間有明顯的相關(guān)性Cr和Ni、Ni和Zn也有比較明顯的相關(guān)性。

表5-2-4交通區(qū)重金屬元素相關(guān)性

在主干道路區(qū)采集的土壤樣品中，Cd、Cr、Cu、Ni、Pb、Zn兩兩之間具有明顯的相關(guān)性，所以可能有相同的污染源。As、Hg與其自身外七種重金屬的相關(guān)性都不明顯，兩者與其他重金屬的來(lái)源可能不同。

表5-2-5公園區(qū)重金屬元素相關(guān)性

在公園綠地區(qū)采集的土壤樣品中，As、Cr、Ni兩兩之間的相關(guān)性比較明顯，可能有相同的污染源。Cu、Pb、Zn兩兩之間的相關(guān)性比較明顯，可能有相同的

14

污染源。

5.2.2、土壤重金屬主成分分析

為了進(jìn)一步明確該城區(qū)不同功能區(qū)土壤中重金屬的污染來(lái)源, 利用Matlab軟件對(duì)上述不同功能區(qū)土壤重金屬的指標(biāo)進(jìn)行主成分分析。

在利用主成分分析進(jìn)行源解析前,采用KMO法和Bartlett法對(duì)原始數(shù)據(jù)集進(jìn)行主成分分析適宜性檢驗(yàn)。經(jīng)計(jì)算, 其KMO值為0.73（大于0. 5的最小值要求），Bartlett檢驗(yàn)統(tǒng)計(jì)的顯著水平小于0. 01,因此可認(rèn)為相關(guān)系數(shù)矩陣與單位陣有顯著差異,以上檢驗(yàn)結(jié)果均表明原始數(shù)據(jù)集適合進(jìn)行主成分分析。

數(shù)據(jù)處理及分析結(jié)果見下表。

表5-2-6生活區(qū)主成分分析 表5-2-7工業(yè)區(qū)主成分分析

生活區(qū)的8種重金屬的全部信息可由3個(gè)主成分表示。即對(duì)前3個(gè)主成分進(jìn)行分析已經(jīng)能夠反映全部數(shù)據(jù)的大部分信息。結(jié)合表中數(shù)據(jù)可以看出:三個(gè)主成分中第一主成分的貢獻(xiàn)率為45.2%，特點(diǎn)表現(xiàn)為因子變量在元素Pb、Cd上有較高的載荷。研究結(jié)果表明汽車尾氣顆粒物中鉛含量較高,在怠速狀態(tài)下向大氣的排放量分別為01571 mg/min。又知Pb和Cd在該區(qū)相關(guān)性明顯，所以大致可認(rèn)為交通污染是生活區(qū)的一個(gè)主要污染原因。第二主成分的貢獻(xiàn)率為14.2%，特點(diǎn)表現(xiàn)為因子變量在元素As、Zn上有較高的載荷，因?yàn)樯罾袝?huì)富含比較多的As和Zn，所以生活區(qū)生活垃圾的堆積可能是重金屬污染的主要原因。第三主成分的貢獻(xiàn)率為13.4%，特點(diǎn)表現(xiàn)為因子變量在元素Ni上有較高的載荷，普遍認(rèn)為，Ni和Zn主要來(lái)源于輪胎的磨損，所以也表明交通污染可能是生活區(qū)的污染來(lái)源之一。

工業(yè)區(qū)的8種重金屬的全部信息可由2個(gè)主成分表示。第一主成分的貢獻(xiàn)率為65.7%，特點(diǎn)表現(xiàn)為因子變量在元素Cr、Cu、Hg、Pb、Zn上有較高的載荷,且這些重金屬在該區(qū)相關(guān)性顯著，可以認(rèn)為來(lái)源相同，工業(yè)污染排放是工業(yè)區(qū)重金屬污染的主要原因。第二主成分的貢獻(xiàn)率為15.8%，特點(diǎn)表現(xiàn)為因子變量在元素As、Hg上有較高的載荷, 這兩種元素通常被視為燃煤和垃圾燃燒的標(biāo)志元素，所以燃煤污染也是工業(yè)區(qū)重金屬污染的一個(gè)主要來(lái)源。

15

表5-2-8山區(qū)主成分分析 表5-2-9主干道路區(qū)主成分分析

山區(qū)的8種重金屬的全部信息可以由3個(gè)主成分表示。第一主成分的貢獻(xiàn)率為38.0%，特點(diǎn)表現(xiàn)為因子變量在元素Zn上有較高的載荷。第二主成分的貢獻(xiàn)率為25.5%，特點(diǎn)表現(xiàn)為因子變量在元素Cd、Pb上有較高的載荷，第三主成分的貢獻(xiàn)率為19.4%，特點(diǎn)表現(xiàn)為因子變量在元素Hg上有較高的載荷。綜合考慮，因?yàn)樯絽^(qū)的污染程度較輕，所以重金屬污染的主要原因可能是山區(qū)富含較多的礦石。

主干道路區(qū)的8種重金屬的全部信息可以由3個(gè)主成分表示。第一主成分的貢獻(xiàn)率為46.9%，特點(diǎn)表現(xiàn)為因子變量在元素Pb、Ni、Cu、Cd上有較高的載荷，所以交通污染中的汽車尾氣排放和輪胎磨損可能是重金屬污染的主要原因。第二主成分中Hg具有較高的載荷，所以也可能受到工業(yè)污染如儀表、電氣工業(yè)等的影響。第三主成分中As具有較高的載荷，所以可能受到工業(yè)上排放As的工廠如化工、冶金業(yè)等的影響。

表5-2-10公園綠地區(qū)主成分

公園綠地區(qū)的8種重金屬的全部信息可以由3個(gè)主成分表示。第一主成分的貢獻(xiàn)率為48.9%，特點(diǎn)表現(xiàn)為因子變量在元素Cd、Cr、Pb、Zn上有較高的載荷，這些重金屬元素主要來(lái)源可以認(rèn)為是交通污染。 第三主成分中Hg具有較高的載荷，有研究發(fā)現(xiàn)由于大量的熒光燈、汞燈等氣體放電燈的廣泛使用,使得綠地表層土壤中汞的含量遠(yuǎn)遠(yuǎn)超過(guò)背景值。所以熒光燈等的大量使用可能是Hg污染的主要原因。

16

5.3 問(wèn)題三模型的建立和求解

5.3.1 污染物傳播特征的分析

重金屬元素傳播途徑主要有水流傳播、土壤滲透、大氣擴(kuò)散等。在城區(qū)我們主要考慮以大氣擴(kuò)散方式傳播。由此，我們通過(guò)擴(kuò)散方程建立模型，從中體現(xiàn)重金屬傳播的主要特性。

5.3.2模型建立和求解

由問(wèn)題分析得，我們主要采用對(duì)流擴(kuò)散方程建立模型。

?c?2c?2c?2c?c?c?c?Ex2?Ey2?Ez2?ux?uy?uz （5?3?1）?t?x?y?z?x?y?z上式

中C為重金屬濃度，Ex，Ey，Ez為在x，y，z三個(gè)方向的擴(kuò)散系數(shù)，ux，uy，uz為在x，y，z三個(gè)方向的對(duì)流系數(shù)。

模型一：基于污染源的（5-3-1）擴(kuò)散模型，建立的逆求污染源模型的基本思想是通過(guò)擴(kuò)散方程離散化，即把連續(xù)的定解區(qū)域用有限個(gè)離散點(diǎn)構(gòu)成的網(wǎng)格來(lái)代替， 這些離散點(diǎn)稱作網(wǎng)格的節(jié)點(diǎn)；把連續(xù)定解區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)來(lái)近似；把原方程和定解條件中的微商用差商來(lái)近似，積分用積分和來(lái)近似，于是原微分方程和定解條件就近似地代之以代數(shù)方程組，即有限差分方程組，解此方程組就可以得到原問(wèn)題在離散點(diǎn)上的近似解。

針對(duì)本問(wèn)題由于污染情況都大致在0~100海拔高度以及忽略對(duì)流影響，以及考慮在x、y方向上的擴(kuò)散系數(shù)相同，得到二維方程：

??2c?2c??c?E??2? 5-3-2 2???t?y???x把上式離散化得到方程：

?c(x??x,y,t)?2c(x,y,t)?c(x??x,y,t)c(x,y??y,t)?2c(x,y,t)?c(x,y??y,t)?c(x,y,t??t)?c(x,y,t)?E*???2?t(?x)(?y)2??

由于污染源點(diǎn)處實(shí)際上具有濃度隨時(shí)間變化相對(duì)較大的情況，問(wèn)題轉(zhuǎn)化為求濃度

c隨時(shí)間t的變化率最大問(wèn)題，進(jìn)一步轉(zhuǎn)化為在等步長(zhǎng)得時(shí)間內(nèi)求濃度c向x,y方

向的擴(kuò)散率。這樣模型可簡(jiǎn)化為只對(duì)相鄰網(wǎng)格求濃度差，通過(guò)比較濃度差的值，來(lái)相應(yīng)確定。

對(duì)模型一的求解，由于本題給出的數(shù)據(jù)具有差異性，這里先利用Matlab對(duì)數(shù)據(jù)插值擬合為后續(xù)的x，y取網(wǎng)格化數(shù)據(jù)提供基礎(chǔ)。這里插值對(duì)本題模型精度造成極大影響。由此我們提出了利用遺傳算法反演擴(kuò)散特征進(jìn)而反求污染源的模型二。

模型二， 在重金屬傳播特征的基礎(chǔ)上引入污染源的確定模型如下：

?C?C?C?C?ux?uy?uz??t?x?y?zEx?C?C?C?E?E??Mi??x?xi???y?yi???z?zi?yz222?x?y?zi?1222q 5-3-3

上式中C為重金屬濃度，Ex，Ey，Ez為在x，y，z三個(gè)方向的擴(kuò)散系數(shù)，ux，uy，

17

uz為在x，y，z三個(gè)方向的對(duì)流系數(shù)。Mi為污染源的排放強(qiáng)度。

?M??x?x???y?y???z?z? 為污染源處對(duì)濃度變化率的影響。進(jìn)而可以通過(guò)

iiiii?1q確定該項(xiàng)的值來(lái)確定污染源的大致位置。

記r??x,y,z?，0<||r||???C?C?2??t?u?r?ECq?r2??Mi??r?ri? 0?||r||?L ,t?0 i?1???C?0,t??0, C?L,t??0 ????C?r,0??0?對(duì)式5-3-4進(jìn)行正向求解如下： 令函數(shù)變換：

??uru2tC?r,t??V?r,t? e??2E-?4E??? 定解問(wèn)題5-3-4可轉(zhuǎn)化為定解問(wèn)題5-3-6

???V?2V??uru2t??2E-4E??q???t?E???r2?e?Mi??r?ri?i?1? ?V?0,t??0, V?L,t??0 ???V?r,0??0那么可得Tn?t?滿足的方程為：

?2?'?n?Tn?t???????L??Tn?t??fn ?Tn?0??0???u2t???n??2?t??B?解得：T?4E???t?L?n??u2?n??2?e?e???? 4E???L???? 18

5-3-4 5-3-5

5-3-6 其中

?in?ri2qB??Mie2EsinLi?1Lux

最后求出方程組5-3-6的解為：

?in?ri2q2E2?2Miesin?n???????ut??t??L?n?rLi?1L???4E?????e?esin2??L 5-3-7 u2?n???????4E?L?uxV?r,t???n?1?為了確定污染源的位置，我們要利用先驗(yàn)數(shù)據(jù)，本文中V|r?ri?V?r,t?為已知。

由5-3-7可知解V?r,t?依賴于ri及Mi,i?1,2,?q，可以把解V?r,t?記為

V?V(r,t,r1,r2,?rq,M1,M2,?Mq)，根據(jù)已知附加條件,源項(xiàng)識(shí)別反問(wèn)題便轉(zhuǎn)化為下面非線性優(yōu)化問(wèn)題：

ri??0,Lmin?V(r,T,r,r,?r,M??12qi?1q1,M2,?Mq)?V(r,T) 5-3-8

?2即使5-3-8最小的ri,Mi,i?1,2,?q就是源項(xiàng)識(shí)別反問(wèn)題的最優(yōu)解。

針對(duì)本問(wèn)題的求解，考慮到對(duì)流對(duì)該模型影響較小，同時(shí)為了簡(jiǎn)化模型方便求解，我們忽略了方程的對(duì)流項(xiàng)，令u=0。此外，由問(wèn)題一的分析得出，重金屬分布主要集中在海拔為0—100m范圍內(nèi)的低海拔區(qū)，在這個(gè)高度范圍內(nèi)我們認(rèn)為在z 方向上沒(méi)有濃度梯度, 即污染物濃度值不變。所以令r=(x,y)，由此我們可以得到式5-3-4、5-3-5、5-3-7的簡(jiǎn)化式如下：

?C?2Cq?E2??Mi??r?ri? ?t?ri?1C?r,t??V?r,t?

V?r,t???n?1?n?ri2q2Msin?n?????i?t??n?rLi?1L??L???1?esin2L ???n???????L?結(jié)合問(wèn)題一得到的各個(gè)重金屬濃度分布圖，先對(duì)位置r歸一化，如||r||=（||r||+1000）./(max(||r||)+1000)。處理通過(guò)Matlab編程，使用遺傳算法，其中參數(shù)取法如見附錄程序，求得各個(gè)重金屬污染源點(diǎn)如下：

（1）As的污染源確定：先對(duì)數(shù)據(jù)表通過(guò)對(duì)As列降序排序，實(shí)際上污染源點(diǎn)必在有污染區(qū)。取其前161組最大濃度點(diǎn)進(jìn)行計(jì)算，再則利用問(wèn)題一的As濃度分布圖，

19

可以預(yù)判污染源個(gè)數(shù)小于等于3個(gè)，以此為基礎(chǔ)提高算法效率。 污染源編號(hào) 1 2 3 污染源位置(x,y) (10856，14727) (6837,3490) (5734，9659) 同理可得： （2）Cd的污染源確定 污染源編號(hào) 1 2 3 4 污染源位置(x,y) (5734，9659) (15748,15728) (205,11228) (9277,16148) （3）Cr的污染源確定 污染源編號(hào) 1 2 污染源位置(x,y) (22965,13535) (8403,1075) （4）Cu的污染源確定 污染源編號(hào) 1 2 污染源位置(x,y) (23359,5325) (10510,15314) （5）Hg的污染源確定 污染源編號(hào) 1 2 3 污染源位置(x,y) (12442，4329) (20261,7586) (10700,8184) （6）Ni的污染源確定 污染源編號(hào) 1 2 3 污染源位置(x,y) (22193，12185) (9277,16148) (13175,12238) （7）Pb的污染源確定 污染源編號(hào) 1 2 3 污染源位置(x,y) (5394，8631) (11902,7709) (10510,15314) （8）Zn的污染源確定 污染源編號(hào) 1 2 3 4 污染源位置(x,y) (7056，8348) (5503,1127) (13244,7056) (14177,6684) 5.3.3模型可行性分析

通過(guò)對(duì)整個(gè)城區(qū)里八種重金屬元素的相關(guān)性分析，可以推出其中有些元素可能具有相同的污染源。這樣可以驗(yàn)證問(wèn)題三的模型是否可行。這里通過(guò)matlab編程得到其相關(guān)性矩陣R為

20

R矩陣橫向排序?yàn)锳s，Cd,，Cr，Cu，Hg，Ni，Pb，Zn；縱向排序?yàn)锳s，Cd,，Cr，Cu，Hg，Ni，Pb，Zn。由R矩陣可知Cr與Ni相關(guān)性系數(shù)為0.7158，Pb與Cd的相關(guān)性系數(shù)為0.6603，理論上可判為有同一污染源點(diǎn)。在問(wèn)題三的求解結(jié)果中Ni的一個(gè)污染點(diǎn)（22193,12185）與Cr的一個(gè)污染源點(diǎn)（205,11228）相近。Pb的一個(gè)污染源點(diǎn)（5394,8631）與Cd的一個(gè)污染源點(diǎn)（5734,9659）相近。符合相關(guān)性的判斷規(guī)律。這證明了模型的可行性與可靠性。

5.4問(wèn)題四的分析和求解

5.4.1優(yōu)缺點(diǎn)分析 優(yōu)點(diǎn)：

本模型參考水流模型、大氣擴(kuò)散模型，比較合理地模擬重金屬在該城區(qū)內(nèi)的擴(kuò)散。通過(guò)對(duì)模型的簡(jiǎn)化，比較簡(jiǎn)單的將模型轉(zhuǎn)化為二維偏微分方程。由于遺傳算法具有極好的全局優(yōu)化性能,減少了陷入局部極值的風(fēng)險(xiǎn),特別對(duì)非線性優(yōu)化問(wèn)題具有很強(qiáng)的適應(yīng)性。傳統(tǒng)的反演方法在源強(qiáng)和濃度的關(guān)系上,只能給出一種虛擬的聯(lián)系,而遺傳算法可以給出較為直接的關(guān)系,可以在很大的程度上,減少對(duì)大量的小源和不確定性污染源調(diào)查的工作量,從而為污染預(yù)報(bào)和污染控制提供了基礎(chǔ)。用遺傳算法來(lái)反演,方法簡(jiǎn)單且通用性強(qiáng)。

缺點(diǎn)：

由于不考慮模型里的對(duì)流項(xiàng)以及海拔上面的影響，導(dǎo)致其自身的確定誤差較大，其次遺傳算法本身就是一個(gè)枚舉算法，算法在運(yùn)行的過(guò)程中具有隨機(jī)性，導(dǎo)致求出的污染點(diǎn)的分布不是完全符合，但具有穩(wěn)定的收斂性，是可以接受的。

5.4.2模型優(yōu)化

在考慮擴(kuò)散方程時(shí)，若能加入對(duì)流項(xiàng)及垂直方向上的擴(kuò)散影響，這樣可以增加擴(kuò)散方程的精度型，進(jìn)而可以增加模型利用遺傳算法的求解精度。其次，如能提供該地區(qū)具體的工業(yè)排放區(qū)及交通要道等的具體密集情況，進(jìn)一步減小污染源位置的范圍和個(gè)數(shù)的初始值，提高遺傳算法得出的污染源的具體位置的準(zhǔn)確性。

21

6、參考文獻(xiàn)

[1] 金昌市城區(qū)土壤重金屬空間分布及潛在生態(tài)危害評(píng)價(jià) 黃璜,南忠仁,胡小娜,劉曉文,李媛,丁海霞 環(huán)境監(jiān)測(cè)管理與技術(shù) 2009年10月 第21卷 第5期

[2] 大氣污染物擴(kuò)散模型應(yīng)用研究 王寒梅, 賈軼然 遼寧師專學(xué)報(bào) 2006年6月 第8卷第2期

[3] 對(duì)流-擴(kuò)散方程源項(xiàng)識(shí)別反問(wèn)題的遺傳算法 閔濤, 周孝德, 張世梅, 馮民權(quán)水動(dòng)力學(xué)研究與進(jìn)展 2004年7月 A 輯第19卷第4 期

[4] 熱傳導(dǎo)方程有限差分法的MATLAB實(shí)現(xiàn) 史策 咸陽(yáng)師范學(xué)院學(xué)報(bào) 2009年7月 第24卷 第4期

[5] 二維穩(wěn)態(tài)對(duì)流—擴(kuò)散方程參數(shù)反演的迭代算法 閔濤，劉相國(guó)，張海燕，艾克鋒 水動(dòng)力學(xué)研究與進(jìn)展 2007年11月 A 輯第22卷第6期

22

附錄

%相關(guān)性分析程序 num=xlsread('a.xls'); L=zeros(319,8); for i=1:8

L(:,i)=num(:,i+1); end

L(:,2)=L(:,2)/1000; L(:,5)=L(:,5)/1000; R1=zeros(8,8); R2=zeros(8,8); R3=zeros(8,8); R4=zeros(8,8); R5=zeros(8,8); for i=1:8 for j=1:8

M=corrcoef(L(1:44,i),L(1:44,j)); R1(i,j)=M(1,2); end end

for i=1:8 for j=1:8

M=corrcoef(L(45:80,i),L(45:80,j)); R2(i,j)=M(1,2); end end

for i=1:8 for j=1:8

M=corrcoef(L(81:146,i),L(81:146,j)); R3(i,j)=M(1,2); end end

for i=1:8 for j=1:8

M=corrcoef(L(147:284,i),L(147:284,j)); R4(i,j)=M(1,2); end end

for i=1:8 for j=1:8

23

M=corrcoef(L(285:319,i),L(285:319,j)); R5(i,j)=M(1,2); end end

[v1 d1]=eig(R1); [v2 d2]=eig(R2); [v3 d3]=eig(R3); [v4 d4]=eig(R4); [v5 d5]=eig(R5); R1 v1 d1 R2 v2 d2 R3 v3 d3 R4 v4 d4 R5 v5 d5

t1=zeros(8,3); t2=zeros(8,3); t3=zeros(8,3); t4=zeros(8,3); t5=zeros(8,3); for i=6:8

t1(:,i-5)=v1(:,i).*sqrt(d1(i,i)); t2(:,i-5)=v2(:,i).*sqrt(d2(i,i)); t3(:,i-5)=v3(:,i).*sqrt(d3(i,i)); t4(:,i-5)=v4(:,i).*sqrt(d4(i,i)); t5(:,i-5)=v5(:,i).*sqrt(d5(i,i)); end t1 t2 t3 t4 t5

24

%遺傳算法V函數(shù)

function [f]=fun(N,V,X) sum=0; q=3;

%q=sum(N); M=10*rand(q,1); x=zeros(q,1); zn=319; T=100; t=1; L=pi; u=0; E=-1; f=0; j=1;

for i=1:siz

if N(:,i)==1 x(j)=X(i);

for n=1:T b=n*pi/L; sum1=0; for k=1:q

m1=M(k).*exp(-u/(2*E).*X(k)).*sin(b.*X(k)); sum1=m1+sum1; end

tmp=2/L*sum1/(b^2)*(exp(0)-exp(-t*(b^2)))*sin(b*x(j)); sum=tmp+sum; end j=j+1; if j==q j=1; end

f=f+(sum-V(i))^2; end

sum=0; end end

25

%遺傳算法主函數(shù)gen.m

um = 100; %種群規(guī)模 Iter = 10; %迭代次數(shù) pinter = 0.9; %交叉概率 pmut = 0.5; %變異概率 %產(chǎn)生初始種群

siz = 161;

for i = 1:Num Temp = rand; Te=3;

% Te = (rand>(1/4)) + (rand>(2/4)) + (rand>(3/4)) + (rand>0) + 1; %5 %Te = (rand>(1/3)) + (rand>(2/3)) + (rand>0) + 1; %4 %Te = (rand>(1/2)) + (rand>0) + 1; %3 %Te = (rand>0) + 1; %2

N(i,:)=(randperm(siz)<=Te); end

igBest = 1; %歷史最佳個(gè)體索引 NBest = N(1,:); %個(gè)體最佳值 F = zeros(1,Num); data=xlsread('a.xls'); V=data(1:siz,2); X=data(1:siz,11);

X=(X+1000)./(max(X)+1000);

%D = xlsread('C:\\Documents and Settings\\Administrator\\桌面\\數(shù)模2011\\最短距離D.xls');

%CAS = xlsread('C:\\Documents and Settings\\Administrator\\桌面\\數(shù)模2011\\A區(qū)發(fā)案率.xls');

%計(jì)算適應(yīng)度 for k = 1:Num

F(k) = fun(N(k,:),V,X); end

it = 1;

while it < Iter

%選擇兩個(gè)個(gè)體交叉，變異 ind = randperm(Num); for k = 1:2:Num %交叉操作

26

n1 = N(ind(k),:); n2 = N(ind(k+1),:);

[n1 n2] = intercross(n1,n2,pinter); %變異操作

n1 = mutate(n1,pmut); n2 = mutate(n2,pmut);

F1 = fun(n1,V,X); F2 = fun(n2,V,X); %較優(yōu)的替換父代

if F1 < F(ind(k)) | F1 < F(ind(k+1)) if F(ind(k)) < F(ind(k+1)) N(ind(k),:) = n1; F(ind(k)) = F1; else

N(ind(k+1),:) = n1; F(ind(k+1)) = F1; end end

if F2 < F(ind(k)) | F2 < F(ind(k+1)) if F(ind(k)) < F(ind(k+1)) N(ind(k),:) = n2; F(ind(k)) = F2; else

N(ind(k+1),:) = n2; F(ind(k+1)) = F2; end

end end

it = it + 1; end

[value ind] = min(F); N(ind,:)

27

%遺傳算法子函數(shù)

function [s1,s2] = intercross(s1,s2,pinter) iTemp = 0; if pinter > rand

ind1 = find(s1); %找到第一個(gè)染色體第k行的1的位置 ind2 = find(s2); %找到第二個(gè)染色體第k行的1的位置

inter = intersect(ind1,ind2); ind1 = setdiff(ind1,inter);

ind2 = setdiff(ind2,inter); %找到兩個(gè)染色體數(shù)字不同的所有位數(shù) if length(ind1) ~= 0 & length(ind2)~=0

%隨機(jī)選擇第一個(gè)染色體的一個(gè)1與另一個(gè)染色體的0互換，然后隨機(jī)選擇第一個(gè)染色體的一個(gè)0（不是剛才的位）跟第二個(gè)染色體的1調(diào)換 iRa = round(rand*(length(ind1)-1)+1); iRa2 = round(rand*(length(ind2)-1)+1); s1(ind1(iRa)) = 1 - s1(ind1(iRa)); s1(ind2(iRa2)) = 1 - s1(ind2(iRa2)); s2(ind2(iRa2)) = 1 - s1(ind2(iRa2)); s2(ind1(iRa)) = 1 - s2(ind1(iRa)); end else

s1 = s1; s2 = s2; end

%變異操作

function s = mutate(s,pmu) PlaneNum = size(s,1); if pmu>rand

ind1 = find(s);

ind2 = find(s ~= 1);

if length(ind1)~=0 & length(ind2) ~= 0 %隨機(jī)找到一個(gè)0跟一個(gè)1的位置，進(jìn)行交換變異 iRa1 = round(rand*(length(ind1)-1))+1; iRa2 = round(rand*(length(ind2)-1))+1; s(ind1(iRa1)) = 1-s(ind1(iRa1)); s(ind2(iRa2)) = 1-s(ind2(iRa2)); end else s = s; end

28