王曉彤，蔣洪迅

(中國人民大學信息學院，北京 100872)

0 引言

東北老工業(yè)區(qū)的空氣污染問題，既影響經(jīng)濟也危害民生，是當?shù)鼐用裆眢w健康的一大隱患.[1-5]在諸多空氣污染物中，細顆粒物PM2.5的影響非常大.細顆粒物PM2.5指環(huán)境空氣中空氣動力學當量直徑小于等于2.5 μm的顆粒物，其在空氣中濃度越大，代表空氣污染程度越高.相較于其他污染物，PM2.5由于其直徑較小，能在大氣中停留較久并隨氣流運動距離較遠，對人類的健康影響巨大[6].由于PM2.5等顆粒物的濃度上升，導致全球每年約210萬人死于心血管和呼吸道疾病(93%)、肺癌(7%)[7].在我國PM2.5已成為最主要空氣污染物之一.因此，PM2.5濃度預測研究，對國計民生具有重要的意義和價值.如果能提供較為準確的預報結(jié)果，不僅可以敦促居民進行提前預防，諸如出行佩戴口罩、減少戶外活動等，還可以通過預測過程建模和數(shù)據(jù)挖掘，探索產(chǎn)生高濃度PM2.5的污染來源和形成條件，為空氣污染防治提供重要的理論依據(jù)和數(shù)據(jù)支持.

目前，針對PM2.5或者其他空氣污染物的預測主要可以劃分為數(shù)據(jù)建模和方法建模兩大類：

數(shù)據(jù)建模是指輸入數(shù)據(jù)的選取和增減，即從數(shù)據(jù)工程或特征工程的角度出發(fā)，研究PM2.5濃度相關的經(jīng)濟、社會、人類行為以及氣象條件的數(shù)據(jù)特征，增減預測模型的輸入數(shù)據(jù)維度和幅度，以探求最優(yōu)的預測效果.現(xiàn)有的空氣質(zhì)量預測，一般來說，都是在選定某種預測模型或者某幾種模型集成的條件，基于現(xiàn)有的歷史空氣質(zhì)量數(shù)據(jù)(如PM2.5，PM10，O3，SO2，CO2等)，再結(jié)合對應時空環(huán)境下的氣象條件數(shù)據(jù)(如溫度、濕度、風速、風向、氣壓、光照時間等)進行機器學習和訓練.

方法建模是指預測方法的調(diào)優(yōu)與改進，即在給定的輸入數(shù)據(jù)不變的情況下，通過對比多個預測模型選擇效果較好的方法，或者通過針對某一給定預測模型的方法改進，或者通過參數(shù)調(diào)優(yōu)獲得更好的預測結(jié)果，抑或集成兩種或兩種以上預測模型以取長補短、提高預測準確度的方法.空氣質(zhì)量預測的方法建模，其主要模型可以歸結(jié)為多元回歸、神經(jīng)網(wǎng)絡和群智能模型等幾大類方法.許多學者采用了多元回歸模型進行了空氣質(zhì)量預測研究[9-11]；自神經(jīng)網(wǎng)絡出現(xiàn)以來，人們開始使用神經(jīng)網(wǎng)絡模型對空氣質(zhì)量進行預測[12-18]；基于進化計算(Evolutionary Algorithm，EA)的群智能模型，也是預測領域一類重要的研究方法.而遺傳規(guī)劃(Genetic Programming，GP)是相對出現(xiàn)較晚的一種EA方法.不少學者將GP及其變體應用在氣象、水文預測等研究領域中，也取得了較好的效果.[19-26]

本文提出一種基于多基因遺傳規(guī)劃(Multi-Gene Genetic Programming，MGGP)的空氣質(zhì)量預測模型，是在經(jīng)典遺傳規(guī)劃基礎上發(fā)展出來的一種群進化智能模型，對于多維時間序列預測問題比較有效且具有更強的魯棒性，特別適合于城市空氣質(zhì)量預測.

1 基于MGGP的空氣質(zhì)量預測模型

遺傳規(guī)劃(Genetic Programming，GP)是遺傳算法(Genetic Algorithm，GA)的一個分支，而MGGP是遺傳規(guī)劃的一種魯棒性變體.為了更好地定義與規(guī)范基于MGGP的空氣質(zhì)量預測模型，首先要限定面向預測領域的經(jīng)典GP模型基本范式和方法框架.

1.1 經(jīng)典遺傳規(guī)劃模型

GP的基本算法框架與GA相同，都是模仿達爾文的進化論，即選擇性地淘汰種群中不適應環(huán)境的個體，并從優(yōu)秀的個體中繁育新一代個體.進化的核心機制是代際更替中的優(yōu)勢個體之間的交叉遺傳以及隨機變異演化，通過一個預先設定的質(zhì)量標準形式評判“種群”解集中每個個體適應度，基于此提供被保留到下一代的概率，經(jīng)過若干次迭代滿足終止條件后，可得到一個最優(yōu)解或近優(yōu)解.一般來說，終止條件既可以是設定的最高進化代數(shù)，也可以是判定的最小誤差等適應度(fitness)標準.

GP跟GA相比不同的是可以生成適用問題模型的啟發(fā)式規(guī)則，或者說是其編碼形式可具解釋性.以最常見的樹形編碼為例，GP的演化運算是對解析樹(parse tree)進行操作，而不是傳統(tǒng)地對比特串(bit string)進行操作.PM2.5濃度預測模型的GP編碼案例見圖1，其編碼形式顯然具有邏輯上的可解釋性，其對應表達式為T(t)/RH(t)-[sin(p(t))+PM(t-1)]，其中T(t)為t時刻的氣溫，RH(t)表示t時刻的相對濕度，p(t)表示t時刻的氣壓，PM(t-1)則是上一時刻的PM2.5濃度.一般來說，解析樹是由一個終止符集(terminal set，TS)和一個函數(shù)集(function set，F(xiàn)S)組成.其中，TS包括函數(shù)的各項參數(shù)和變量(常數(shù)、邏輯常量、變量等)；FS包括基本算術運算、標準編程操作、標準數(shù)學函數(shù)、邏輯函數(shù)或其他數(shù)學函數(shù).

已有的研究表明，PM2.5是由空氣中原有污染物在一定的氣象條件下，經(jīng)過一系列未知的物理化學過程形成的氣溶膠性質(zhì)懸浮顆粒物.目前已知影響因素中，PM2.5濃度與SO2和NO2等污染物排放及濕度、日照、風力、風向等氣象條件有關，甚至與人類活動、地形環(huán)境等諸多因素都相關.由于遺傳規(guī)劃GP不僅在復雜非線性回歸空間上具有高成功率，而且能夠生成適用問題模型的啟發(fā)式規(guī)則表達式，從而為一些潛在的過程提供一些可能的解釋，所以GP在空氣質(zhì)量預測領域中特別是PM2.5濃度預測中具有明顯的優(yōu)勢.

1.2 MGGP模型

MGGP是傳統(tǒng)遺傳規(guī)劃GP的一種變體，即通過線性組合低深度的GP樹來提高傳統(tǒng)GP的適應性和魯棒性.[27]由于使用較小的解析樹，MGGP有望提供比傳統(tǒng)GP更簡單的模型[28]. 在MGGP中，預測變量通過多基因程序中的每個基因(即解析樹)的加權輸出加上偏移項來計算. 例如，MGGP個體使用包括t時刻的氣溫、相對濕度、氣壓，上一時刻PM2.5濃度，通過其輸入數(shù)據(jù)的2個基因來預測t時刻PM2.5濃度(見圖2)，在數(shù)學上，這個模型可以寫成

(1)

其中d0為偏移項，d1，d2為回歸系數(shù)(即每個基因的權重).通常，系數(shù)由每個MGGP個體的普通最小二乘法確定[29].因此，MGGP采用經(jīng)典線性回歸方法來捕獲非線性行為，而不需要預先指定非線性結(jié)構.

在MGGP中，每個基因都是一個簡單的GP解析樹，不會任何隱含或明確地參考同一染色體中的其他基因.每個初始染色體可能包含一個或幾個標準的低深度GP解析樹，染色體中基因的最大數(shù)量和基因的深度可自行定義.然后，使用標準的子樹交叉和變異轉(zhuǎn)換以及直接復制來產(chǎn)生后代.在MGGP中，除了經(jīng)典的GP子樹交叉算子之外，還可以使用被稱為“兩點高級交叉”的樹交叉算子來交換染色體之間的基因.可以設置每個GP運算符的相對概率.

圖1 GP的樹形遺傳編碼案例 圖2 MGGP個體舉例

1.3 基于MGGP的PM2.5濃度預測模型

對于典型回歸方法，其模型結(jié)構和系數(shù)通常在初始階段就需確定，而在初始階段如何合理有效確定這些因素卻非常困難.而對于MGGP，僅有TS和FS需要初始化，學習方法會自主地找到模型的最佳形式和參數(shù).而且MGGP還有一個優(yōu)勢，它能夠自動從輸入變量中尋找對模型有利的變量加以使用，并忽略對模型建立沒有幫助的變量，有效地減少模型的維度，同時增強模型的可理解性.

1.3.1 變量集TS的設定

考慮因素相關性和數(shù)據(jù)可獲得性，本文選取沈陽地區(qū)11個空氣質(zhì)量監(jiān)測站的PM2.5濃度歷史數(shù)據(jù)、相應時空環(huán)境下的氣象條件數(shù)據(jù)作為MGGP模型輸入數(shù)據(jù).各類數(shù)據(jù)的類型、符號和單位，如表1所示.

表1 氣象學數(shù)據(jù)及其符號表示

RH即單位體積空氣中含有的水汽密度與同溫度下飽和的水汽密度的百分比.相對濕度越大，則空氣中所含水蒸氣越多.由于相對濕度為百分數(shù)，故本文用RH(t)表示t時刻的相對濕度的100倍；風向(Wind Direction，WD)是指風吹來的風向，一般以角度表示，正北方向為0°(或360°)，正東方向為90°，正南方向為180°，正西方向為270°，其余風向均由此計算得出，以WD(t)表示t時刻的風向.顯而易見的是前一時段的PM2.5濃度與下一時段的PM2.5濃度有著較強相關性，因此，將前一時段PM2.5濃度也作為輸入數(shù)據(jù)之一；用PM(t-1)表示t-1時刻的PM2.5濃度.在后續(xù)實驗中，實現(xiàn)提前1 d(24 h)、提前2 d(48 h)、提前3 d(72 h)預測當前的PM2.5濃度，表2中提到的t-1時刻在不同實驗中分別代表24，48及72 h前的PM2.5濃度.

1.3.2 FS的設定

在函數(shù)集中，除了基本的數(shù)學四則預算以及IF-THEN條件函數(shù)以外，新模型還引入了一元運算符和三元運算符等較為復雜的計算函數(shù).如SQRT表示平方，EXP表示以自然常數(shù)e為底的指數(shù)函數(shù)，ADD3表示三元加法，即3個變量或常數(shù)連續(xù)相加，MULT3表示三元乘法，即3個變量或常數(shù)連續(xù)相乘.

1.3.3 基因數(shù)的設定

影響MGGP算法性能的相關事項中，一個最主要因素是基因數(shù)目.為摸清基因數(shù)目對模型預測準確性和算法時間復雜度的影響，本文針對基因數(shù)進行了參數(shù)嘗試的預實驗(見圖3).結(jié)果表明，在提前24 h的PM2.5濃度預測中，基因數(shù)目小于5時，均方根誤差RMSE隨基因數(shù)增加而不斷降低；基因數(shù)目大于等于5的情況下，RMSE變化不再是基因數(shù)目的單調(diào)減函數(shù)，而是偶有反彈波動；隨著基因數(shù)目的增加，RMSE下降速度逐漸變緩.類似地，在提前48和72 h的濃度預測中，RMSE分別在基因數(shù)為6和4時開始上下波動，不再單調(diào)遞減，且下降速度變緩.

圖4描述了不同基因數(shù)MGGP模型所需計算時長的變化.顯然，隨基因數(shù)目不斷增加，MGGP所用時間整體不斷延長.其中，提前24 h的預實驗所用時間在基因數(shù)為3和8時均有顯著的增加；提前48和72 h的預實驗時長則分別在基因數(shù)為6和7時開始激增；當基因數(shù)為10時，3個時段的PM2.5濃度預測所需時長已近似達到基因數(shù)目為6時實驗所用時間的2倍、基因數(shù)目為2時實驗所用時間的4倍.

圖3 不同基因數(shù)目對模型預測準確度的影響 圖4 不同基因數(shù)目計算時長

綜上所述，同時考慮模型準確性和算法復雜性，基因數(shù)為6時RMSE相較于基因數(shù)為2時已有顯著下降，而后下降變緩或反彈，且基因數(shù)為6時實驗用時顯著低于基因數(shù)為7～10的情況，故本文選取6作為模型基因數(shù)設定.

1.3.4 MGGP模型其他參數(shù)的設定

本文提出的MGGP空氣質(zhì)量預測模型的其他各項參數(shù)設置見表2.其中：競爭規(guī)模是指在每一代種群中選擇參與競爭的個體數(shù)；競爭規(guī)模越大則適應度較弱的個體被選擇參與遺傳的可能性越小，本文將競爭規(guī)模設為15，即在每代種群中選擇15個個體進行競爭；精英比例表示每一代中可不進行改變而直接進入下一代的最優(yōu)個體比例，本文將精英比例設為0.1，則在每代的1 000個個體中，適應度最優(yōu)的前100個個體可直接進入下一代，而不必進行變換；遺傳操作比例是指參與遺傳的個體使用不同遺傳操作的比例(亦即概率)，本文將變異比例設為14%，交叉比例設為84%，復制比例設為2%.本文使用均方根誤差RMSE作為個體適應度衡量標準，RMSE越低，表示預測回歸誤差越小，適應度越高，本文設置當RMSE降至0.003及以下時，停止進化，否則一直進化至500代.

表2 多基因遺傳規(guī)劃參數(shù)設置

在代碼實現(xiàn)方面采用GPTIPS2.0工具包作為MGGP的求解平臺.GPTIPS2.0是一個基于MATLAB的可完成多基因遺傳規(guī)劃的開源工具包，具有很好的通用性和可擴展性，且不需要其他MATLAB工具包來輔助完成MGGP的實現(xiàn).[28]

2 實證研究

2.1 實驗設置

本文選取沈陽地區(qū)11個空氣質(zhì)量監(jiān)測站，自2016年1月1日0：00到2016年3月26日8：00每小時的各項氣象數(shù)據(jù)及對應PM2.5濃度共22 528組，刪去其中信息殘缺非常嚴重的107組數(shù)據(jù)，使用剩余22 421組數(shù)據(jù)進行實驗.從 22 421組數(shù)據(jù)中隨機抽取20 179組作為訓練集，剩余2 242組作為測試集.

在對比實驗中，本文將選取BP神經(jīng)網(wǎng)絡和傳統(tǒng)GP遺傳規(guī)劃模型作為對照組模型，來衡量MGGP在PM2.5濃度預測中的實際效果.神經(jīng)網(wǎng)絡在PM2.5濃度預測中的表現(xiàn)普遍優(yōu)于多元線性回歸，且BP神經(jīng)網(wǎng)絡作為一種多層前饋神經(jīng)網(wǎng)絡，近年來因其較好的穩(wěn)定性和準確性而較多地應用于預測領域，并被不斷改良使用.[31-34]本文分別使用MGGP和BP神經(jīng)網(wǎng)絡對PM2.5濃度進行全站點和分站點的提前24，48和72 h的預測，并對預測誤差RMSE進行比較分析.另外，為體現(xiàn)MGGP相較于傳統(tǒng)GP的優(yōu)勢，也使用了GP的PM2.5濃度預測作為對照實驗，其參數(shù)設置除不使用多個基因外，其他均與MGGP相同.

2.2 MGGP實驗過程

針對11個站點的全部數(shù)據(jù)，在提前24 h的PM2.5濃度預測實驗中，種群最優(yōu)適應度及平均適應度變化見圖5.初始種群中，最優(yōu)適應度為40.149 3，平均適應度為49.601 1，隨著種群不斷進化，適應度(即均方根誤差RMSE)不斷下降，由于初始種群是有程序隨機生成，誤差較大，因此初始幾代進化較快，誤差迅速減少，隨后進化速度漸緩，進化500代后得到最優(yōu)適應度為38.188.同時，MGGP在進化時優(yōu)先滿足每代種群的最優(yōu)適應度不斷降低，因此在整個進化過程中，各代種群的平均適應度并非是嚴格的單調(diào)遞減函數(shù)，且收斂相對于最優(yōu)適應度稍緩，但整體仍呈下降趨勢，說明種群在進化過程中不斷優(yōu)化.計算公式為：

-0.000 454 7(x4-1.0x5+x6+cos(x2)+10.69)2-186.6；

(2)

(3)

(4)

0.004 151x2(x2+x6+x3+x4+10.87)+0.002 075x5(x2+10.94)(2.0x2+x4)；

(5)

-0.000 526 3(x2+10.94)(x3+2.0x5)(x3+2.0x5+2.0cos(x2))；

(6)

(7)

最優(yōu)個體計算公式為：

(8)

表3自變量及其含義對應情況

自變量變量含義 x1p(t)x2T(t)x3RH(t)x4WD(t)x5WS(t)x6PM(t-1)

公式(2)—(7)分別為最優(yōu)個體的6個基因及偏差項，其中偏差項與第1個基因同在公式(2)中體現(xiàn)，6個基因組成的最優(yōu)個體用公式(8)表示，其中自變量與具體輸入變量含義對應情況見表3.

MGGP提前48 h的PM2.5濃度預測實驗中，種群最優(yōu)適應度及平均適應度變化曲線見圖6，初始最優(yōu)適應度為41.282 5，平均適應度為51.001 2，進化過程中適應度變化趨勢與提前24 h的PM2.5濃度預測實驗中變化趨勢類似，進化500代后最優(yōu)個體的適應度為38.925 2，略高于提前24 h的PM2.5預測情況.最優(yōu)個體的6個基因及偏差項計算公式為：

0.000 522 2-x2(x1-1.0x2)(x4-1.0x6)-18.95；

(9)

(10)

(11)

11.79sin(sin(cos(x5)))；

(12)

0.545 9x2(x4-1.0x6)；

(13)

0.041 37(x3-x5)2+0.041 37x1-0.082 74x4+0.082 74x6+0.279 5.

(14)

最優(yōu)個體計算公式為：

(15)

圖5 MGGP提前24 h PM2.5濃度預測情況 圖6 MGGP提前48 h PM2.5濃度預測情況

各自變量及其對應含義見表3.提前72 h的PM2.5濃度預測中，種群最優(yōu)適應度及平均適應度變化情況見圖7，初始最優(yōu)適應度為41.030 2，平均適應度為50.659 9，進化過程中適應度變化趨勢與前兩項實驗中變化趨勢類似，進化500代后最優(yōu)個體的適應度為38.829 8.最優(yōu)個體的6個基因及偏差項計算公式為：

(16)

(17)

(18)

0.000 289 8(x2+x5)(2.0x2+x3)(2.0x2+x4-x6-7.33)；

(19)

(20)

(21)

最優(yōu)個體計算公式為：

(22)

2.3 對比實驗結(jié)果

本文使用相同數(shù)據(jù)對BP神經(jīng)網(wǎng)絡和傳統(tǒng)遺傳規(guī)劃GP進行預測訓練，得到提前24，48及72 h預測PM2.5濃度最小誤差與對應MGGP最小誤差進行對比(見表4).從表4可以看出，在PM2.5 濃度預測效果上，MGGP與BP神經(jīng)網(wǎng)絡的預測效果受預測時間間隔的影響較小，在提前24，48或72 h的預測中，均能保證預測誤差的相對穩(wěn)定，相比之下，傳統(tǒng)GP則受預測時間間隔影響更大，預測誤差隨時間間隔增大而不斷上升.MGGP與BP神經(jīng)網(wǎng)絡及GP相比體現(xiàn)出了更好的性能，在提前24，48，72 h的濃度預測中產(chǎn)生的誤差，MGGP比BP神經(jīng)網(wǎng)絡降低了4%～10%；而傳統(tǒng)GP在提前24 h的預測中，誤差低于BP神經(jīng)網(wǎng)絡，但在提前48和72 h的預測中，誤差明顯提高，并高于BP神經(jīng)網(wǎng)絡.此外，MGGP在不同時段PM2.5濃度預測中所產(chǎn)生的誤差的標準差比BP神經(jīng)網(wǎng)絡與傳統(tǒng)GP降低了13%～63%，體現(xiàn)了MGGP在不同時段PM2.5濃度預測中的魯棒性優(yōu)于BP神經(jīng)網(wǎng)絡及傳統(tǒng)GP.

表4 MGGP與BP神經(jīng)網(wǎng)絡預測情況對比

MGGP與GP種群進化收斂情況見圖8.GP可在100代以內(nèi)完成收斂，100代之后RMSE下降緩慢，而MGGP從200代后RMSE才開始下降變緩，這是由于MGGP中每個個體由多個基因(解析樹)組成種群復雜度高，因此相比于GP收斂稍慢.但因其基因的多樣性，雖然GP與MGGP的初始種群均由系統(tǒng)隨機產(chǎn)生，MGGP初始的最優(yōu)適應度卻均低于GP，收斂過后3 h段的最終預測準確性均明顯高于GP，體現(xiàn)了多基因的種群優(yōu)勢.同時，隨著預測時間間隔的增加，二者性能差距不斷擴大，MGGP在不同時段的預測中，誤差變化較小，尤其在提前48和72 h的預測中，誤差相差更小，且提前72 h的預測誤差略低于提前48 h的預測誤差，說明當預測時間間隔較大時，MGGP幾乎不受時間間隔變化的影響；而GP提前48和72 h的預測誤差明顯高于提前24 h的誤差，且提前72與48 h的PM2.5濃度預測誤差也相差較大，說明使用MGGP預測PM2.5濃度的魯棒性也顯著優(yōu)于GP.

圖7 提前72 h PM2.5濃度預測情況 圖8 MGGP與GP種群進化情況

針對區(qū)分站點的PM2.5濃度預測，各個站點數(shù)據(jù)使用不同預測模型分別提前24，48，72 h的PM2.5 濃度預測結(jié)果見表5，其中黑體部分數(shù)據(jù)表示該站點在該預測時段中不同模型預測結(jié)果相比的最佳表現(xiàn).從表5可以看出，在11個站點3個時段預測的33個預測情形中，MGGP在其中32個情形下的預測誤差低于GP及BP神經(jīng)網(wǎng)絡，且穩(wěn)定性明顯優(yōu)于GP和BP神經(jīng)網(wǎng)絡.同時，在分站點預測的實驗中，提前24，48，72 h的預測誤差上，GP均大多高于BP神經(jīng)網(wǎng)絡，與全站點中的情況不同；這可能是由于分站點后訓練數(shù)據(jù)顯著減少，導致GP生成的預測模型表現(xiàn)明顯下降，這也從側(cè)面體現(xiàn)了MGGP比GP的穩(wěn)定.此外，BP神經(jīng)網(wǎng)絡在分站點的預測實驗中偶爾會出現(xiàn)較高的誤差，如渾南東路站點提前48 h、京沈路提前72 h及新秀街提前72 h的PM2.5預測誤差，RMSE達到了80～90，同樣從側(cè)面反映出MGGP較高的穩(wěn)定性.

表5 不同站點中不同模型的預測結(jié)果對比

3 結(jié)論

本文提取氣象條件數(shù)據(jù)及PM2.5濃度歷史數(shù)據(jù)，采用多基因遺傳規(guī)劃MGGP對PM2.5濃度進行了全站點及分站點分別提前24，48和72 h的預測，并與BP神經(jīng)網(wǎng)絡及傳統(tǒng)GP的預測水平進行了比較.實驗結(jié)果表明，MGGP在PM2.5預測上具有較好性能，且預測誤差不受時段間隔長短的影響，在不同時段的預測中誤差近似，不會因時間間隔的增加而顯著降低預測的準確率.相較于BP神經(jīng)網(wǎng)絡，MGGP在性能上具有顯著優(yōu)勢，在提前24，48和72 h的預測中均有更小的均方根誤差，且在不同站點不同預測時段的誤差也絕大多數(shù)低于BP神經(jīng)網(wǎng)絡及GP，反映了MGGP對于預測時長變化及空間變化具有更好的魯棒性.但由于基因數(shù)目更多，MGGP比傳統(tǒng)GP收斂速度略低.