MATLAB程序在對污染物主成分分析及預測中的應用

鄭媛媛

（深圳市環境科學研究院，廣東 深圳518006）

1 引言

在近幾年中，城市進行了大量的舊房拆遷及擴建工程，使降塵量增大，可吸入顆粒物也不斷逐步向上發展，造成了濃度的增加。所以，利用2001～2006年的主要污染物濃度和MATLAB技術來分析和預測大氣環境主要污染物的變化趨勢，基礎數據見表1。

表1 2001～2006年空氣主要污染物濃度

2 主成份分析的MATLAB的實現

2.1 使用協方差矩陣進行主成分分析

［PC，LATENT，EXPLAINED］＝pcacov（X）

式中，X為原始指標矩陣的協方關矩陣；PC為主成分（按列放置）；LATENT為協方差矩陣的特征值；EXPLAINED為每個主成分解釋的方差量占總方差量的百分比。

2.2 直接利用原始為進行主成分分析

以下式中，X為原始數據矩陣（每一行代表一個樣本，每一列代表一個指標）；返回SCORE為Z分數；返回TSQUARE為Hotelling的T2統計量；Z分數是通過將原始數據轉換到主成分空間中得到的數據，LATENT向量的值為SCORE的列數據的方差；Hotelling的T2統計量為來自數據集合中心的每一個觀測量的多變量距離的度量。

具體分析的編程如下：

%c9fun9＿1

clear all；clc

%將原始數據樣本矩陣賦值給X；

covX＝cov（X）；%得到原始數據矩陣的協方差矩陣

［PC，LATENT，EXPLAINED］ ＝ pcacov（covX）；%主成分分析

disp（＇綜合因子＇）；PC

disp（＇特征值＇）；LATENT

disp（＇某一綜合因子解釋方差占總方差的百分比%＇）；EXPLAINED

［n，m］＝size（X）；%n為樣本數；m為指標數

Y＝X×PC%n個樣本的m個綜合因子（PC）得分

Y1＝X（：，：）×PC（：，1）%n個樣本的第1個綜合因子（PC）得分

Y2＝X（：，：）×PC（：，2）%n個樣本的第2個綜合因子（PC）得分

plot（Y1，Y2，＇.＇）；xlabel（＇第一主成分＇）；ylabel（＇第二主成分＇）；title（＇樣本分類圖＇）；

for ii＝1∶1∶n；text（Y1（ii），Y2（ii），［＇s＇，num2str（ii）］）；end

綜合因子

特征值

某一綜合因子解釋方差占總方差的百分比%

從主成分分析圖1中可看出：S2（二氧化氮）、S7（硫酸鹽化速率）、S5（空氣污染指數）、S8（pH值）占主要成份，是這幾年對大同市環境污染的主要因素。那么未來是怎樣的變化趨勢呢？為此，對這八項污染因子作了以下進一步的預測分析。

圖1 空氣主要污染物樣本分類

3 用灰色系統預測模型及其在MATLAB的實現

（1）對二氧化硫（mg／m3）二氧化氮（mg／m3）總懸浮物（mg／m3）降塵濃度（t／km2·30d）的預測分析：編制的MATLAB的程序如下。

disp（＇估計出的系統預測模型系統－微分方程組系數＇）

AA＝A

UU＝U（∶）

%調用微分方程組的求解命令求解微分方程組t0＝2001；tf＝2015；

y0＝［15.3 4.5 52 39.52］；

［tt，xx］＝ode15s（＠c4fun44，［t0，tf］，y0，［］）；

plot（tt，xx（∶，1），＇y－∧＇，tt，xx（∶，2），＇b－o＇，tt，xx（∶，3），＇g－＊＇，tt，xx（∶，4），＇r－p＇）；

legend（＇二氧化硫－預測＇，＇二氧化氮－預測＇，＇總懸浮物－預測＇，＇降塵濃度－預測＇）；

xlabel（＇年度＇）；ylabel（＇預測值＇）；hold on；

function dxt＝c4fun44（t，x）

global AA UU；

dxt＝AA＊［x（1）；x（2）；x（3）；x（4）］＋UU；

估計出的系統預測模型系統－微分方程組系數

從圖2可看出降塵濃度有增長趨勢。

圖2 空氣主要污染物濃度預測圖a

（2）對污染指數、可吸入顆粒物、硫酸鹽化速率pH值的預測分析：

disp（＇估計出的系統預測模型系統－微分方程組系數＇）

AA＝A

UU＝U（∶）

%調用微分方程組的求解命令求解微分方程組

t0＝2001；tf＝2015；

y0＝［5.7 1.78 0.899 6.67］；

［tt，xx］＝ode15s（＠c6fun64，［t0，tf］，y0，［］）；

plot（tt，xx（∶，1），＇y－∧＇，tt，xx（∶，2），＇b－o＇，tt，xx（∶，3），＇g－＊＇，tt，xx（∶，4），＇r－p＇）；

legend（＇污染指數－預測，可吸入顆粒物－預測，硫酸鹽化速率－預測，pH值－預測）；

xlabel（＇年度＇）；ylabel（＇預測值＇）；hold on；

function dxt＝c6fun64（t，x）

global AA UU；

dxt＝AA＊［x（1）；x（2）；x（3）；x（4）］＋UU；

c6fun6＿4

估計出的系統預測模型系統－微分方程組系數

從圖3可看出鹽化速率和可吸入顆粒物在逐步向上發展。

圖3 空氣主要污染物濃度預測圖b

4 結語

根據圖1主成分的表明，S2（二氧化氮）、S7（硫酸鹽化速率）占主要污染物，預測圖2、圖3發現降塵濃度、鹽化速率和可吸入顆粒物（PM10）的有增長趨勢。但預測中二氧化氮在今后的變化趨勢中并不明顯，不作為分析的對象。

可吸入顆粒物（PM10）是最近環境各類媒體上，最熱門的話題，PM10污染的形成不外乎兩個途徑。即污染源的貢獻及自然條件，這些因素相互影響共同作用造成PM10污染。這幾年PM10值的上升，與某市車輛的增加和城市的拆遷、擴建改造工程的加大也有一定的關系，根據有關資料表明，該市的2006年機動車數量為2.8063萬輛；2007年機動車數量為3.1049萬輛；2008年是3.483萬輛。預測2012年車輛的可增達30多萬輛。在近兩年中，城市進行了大量的舊房拆遷及擴建工程，使降塵量增大，可吸入顆粒物也不斷逐步向上發展，造成了濃度的增加。

城市的主要污染物是長期以來對城市居民健康造成直接的危害性，為更好地反映其污染變化趨勢、加強污染防治工作和預防嚴重污染事件的發生，研究污染預測方法，這樣的預測將給政府提供了科學有力的治理依據。可更好地開展污染預報及治理工作。為此，借助于MATLAB極強的非線性處理能力，來解決大氣環境主要污染物自身復雜多變的問題，利用MATLAB技術應用到空氣污染的主份因子分析及預測領域中，是很好的治理手段。

［1］ 宋新山，鄧 偉，張 琳.MATLAB在環境科學中的應用［M］.北京：化學工業出版社，2008.

［2］ 吳 萍，陳作平，曾萬明.汽車污染物排放的危害有對策［J］.科技視野，2008（9）：53.

［3］ 劉登國，居 力.機動車排氣污染新標準頒布標準對管理模式變革的探討［J］.環境監測管理與技術，2008（2）：5～6.