季 冕，程 龍

(安徽省環(huán)境監(jiān)測中心站，安徽 合肥 230071)

臭氧(O3)是大氣中常見的微量氣體，主要分布在10 km至50 km的平流大氣層中.高濃度的臭氧會刺激人體組織黏膜，對人體造成傷害.隨著工業(yè)化發(fā)展進程的加快，臭氧污染已經(jīng)成為一個嚴峻的社會問題，甚至在某些城市取代PM2.5成為空氣污染的“罪魁禍首”[1].近地大氣中臭氧的濃度變化有很顯著的季節(jié)規(guī)律，很大程度上也有人為因素的影響，這使每日臭氧濃度變化具有很大的不確定性.

空氣質(zhì)量模型的研究發(fā)展經(jīng)歷了3個階段[2]，第一代的箱式模型、高斯煙團等，第二代考慮復雜擴散過程的歐拉網(wǎng)格模型，到多區(qū)域多尺度的綜合空氣質(zhì)量模型.目前，對臭氧等污染物的濃度預報的第三代模型，主要有兩種[3]：一種是以大氣動力學為基礎(chǔ)，建立關(guān)于大氣污染物濃度稀釋擴散的數(shù)值模型，通過計算來模擬和預測大氣污染物的動態(tài)分布，亦即數(shù)值預報；另一種是基于統(tǒng)計學方法，建立污染物濃度與氣象參數(shù)間的統(tǒng)計預報模型，用來預測大氣污染濃度，亦即統(tǒng)計預報.

目前，主流的空氣預報模式主要有：美國環(huán)境技術(shù)公司(ENVIRON)提出的CAMx[4]、美國國家環(huán)境保護局(USEPA)[3]開發(fā)的第三代空氣質(zhì)量預報和評估系統(tǒng)(Models-3)中的CMAQ(community multiscale air quality model)、中科院大氣物理研究所開發(fā)的嵌套網(wǎng)格空氣質(zhì)量預報模式NAQP(nested air quality prediction)[5]、基于美國環(huán)境預測中心(NCEP)和美國國家大氣研究中心(NCAR)等科研機構(gòu)開發(fā)的WRF(weather research forecast)的WRFC模式[3]等.

1 研究方法和數(shù)據(jù)來源

使用某種單一的模式獨立進行氣象數(shù)據(jù)預測，通常會帶來很大的不確定性，預測偏差和方差往往比較大[6].因為每種方法有其特殊的適應(yīng)性，為方便結(jié)合各模式的優(yōu)點，很多方法將這些模式組合起來，構(gòu)成一個新的預測模式，即集合預報模式.大氣科學界比較著名的最優(yōu)化集合預測方法(operational consensus forecasts, 簡稱OCF)是Woodcock和Engel提出的一種自動化的集合預測系統(tǒng)[7].通常選取一定時間作為預測的滑動窗口，對于一些基礎(chǔ)模式的預測數(shù)據(jù)進行整合升級.

筆者以多個氣象模式預測數(shù)據(jù)和真實觀測數(shù)據(jù)為樣本，利用已有預報模型的某種線性組合或凸組合來構(gòu)成一種全新的臭氧預報模型，并對其性能進行量化評測.總體來講，主要工作如下：首先實現(xiàn)基于Boosting的集合預測方法嶺回歸(ridge regression，簡稱RR)算法，針對算法中參數(shù)的特點進行分析，并給出參數(shù)的調(diào)整方案.其次進行算法性能的測試和評價.根據(jù)我國站點觀測值和模式預報值，對OCF算法和RR算法分別給出相應(yīng)的集合預報值，并通過均方根誤差(RMSE)和時間序列進行性能評估.

原始數(shù)據(jù)來自中國環(huán)境監(jiān)測總站的環(huán)境空氣質(zhì)量數(shù)值預報模式系統(tǒng)，原始數(shù)據(jù)集包含2015年9月1日至2017年2月19日的各個模式的預報結(jié)果和相應(yīng)的觀測值資料.

2 Ridge Regression算法介紹及實現(xiàn)策略

機器學習方法是人工智能領(lǐng)域研究的重要內(nèi)容，集成機器學習的相關(guān)算法往往可以根據(jù)有限的觀測數(shù)據(jù)，繞過難以通過數(shù)值方法直接得出的公式描述，得到所需的結(jié)果數(shù)據(jù)[8-9].筆者介紹的機器學習集合預測方法，可以很好地集成這些模式的優(yōu)勢所在，從而給出一個更加有效的預測結(jié)果.基于PAC可學習性原理中的Boosting算法[10-11]，該文給出了一個較好的思路.

機器學習旨在實現(xiàn)自動決策或預測.通常，需要基于初步估計步驟做出良好決策.然而，并不是所有算法都依賴于估計統(tǒng)計參數(shù)，筆者應(yīng)用的序列集成技術(shù)即不采用估計的算法.

2.1 概 述

2.1.1 符號定義

2.1.2 序列集成技術(shù)

集成學習方法中，組合中每個模式獲得的權(quán)重與模式過去預測值與觀測值的差距有關(guān)系，也就是說理論上預測效果越好的模式將獲得更高的權(quán)重.序列集成技術(shù)類似一個黑盒子，根據(jù)PAC可學習性理論，集成后的預測結(jié)果會優(yōu)于單個預測模式(基學習器).單個的預測模式可以來自統(tǒng)計建模、數(shù)值模式或者決策樹、神經(jīng)網(wǎng)絡(luò)等復雜的方法.所以，該方法實際是工作在元模型級別的.

(1)

2.2 算法詳細說明

RR算法即嶺回歸算法[8]，通常設(shè)置1個參數(shù)λ≥0，并初始化向量u1=(0,…,0)，對于t≥2有

(2)

RR算法流程如下：

fort=1,2,…,T:

輸出：uT=(uT[1],…,uT[N])

2.3 算法實現(xiàn)策略

RR算法實現(xiàn)的關(guān)鍵是通過迭代來不斷更新模型權(quán)重，算法的核心代碼如圖1所示.在確定參數(shù)λ的值時，通?？梢猿跏蓟粋€遍歷范圍，然后以一定的步長進行遍歷，找到給定精確度下最優(yōu)的參數(shù)取值范圍.

算法核心模塊的輸入可以是4個模式矩陣和1個觀測矩陣，其格式可以為s×t(站點數(shù)、日期數(shù))、罰分因子λ.算法輸出為計算出的集成模式數(shù)值，其規(guī)模也為s×t，或者進一步根據(jù)實際觀測矩陣(向量)給出該集成模式的均方根誤差值(RMSE).

算法流程如下：

輸入：P1、P2、P3、P4：4個模式矩陣，OB:觀測矩陣，s：矩陣站點數(shù)，t：日期數(shù)，u：可初始化為0，保存加權(quán)向量，R：保存計算結(jié)果，A：保存每一步的中間變量Ai，AP:保存每一步Ai的廣義逆矩陣，λ：罰分因子

fori=1,…,t:

v=u[i];#保存當天加權(quán)向量

forj=1,…,s

p=[P1[j,i];P2[j,i];P3[j,i];P4[j,i]]; #保存模型在特定t和s下的預測值

R[j,i]=v*p;

Re=R[j,i]-OB[j,i];

e_sum=e_sum+Re*p;

pd=p*p的逆;

pd_sum=pd_sum+pd;

end for

A[i+1]=A[i]+pd_sum; #用每天累加后的pd_sum更新A[i+1]

AP[i+1]=A[i+1]的廣義逆矩陣；

u[i+1]=u[i]-(AP[i+1]*e_sum)的逆陣；

end for

OBA=t-30至t的OB值；

RA=t-30至t計算出的RR值；

輸出：該參數(shù)λ下的RR矩陣R和其RMSE值res.

外層可以通過一個腳本方法，通過步長和區(qū)間的設(shè)置，尋找到合適的罰分因子λ.如果最后輸出的是RMSE值，可以將每次算出的RMSE值存儲到1個數(shù)組中，然后求得其中最小的RMSE罰分因子的λ值.RR算法的性能具有理論上的保證，算法性能的評測標準在下一部分詳細說明.

3 性能評測與分析

3.1 評測指標與度量標準

均方根誤差(root mean square error，簡稱RMSE)亦稱為標準誤差，其定義為：觀測值與真實值差值的平方除以觀測次數(shù)的平方根.它對一組測量中特大值和特小值比較敏感，故常用來反映精度.其表達式如下

(3)

其中：n為觀測次數(shù)，xi為第i次的模式值，yi為第i次的真實值.RMSE值越小，說明模式的精度越高，與真實值貼近的越好.

(4)

其中：t0是評估開始時的第1個時間，T是評估結(jié)束的時間.

針對某個觀測數(shù)據(jù)集和模式預測集，有以下幾個度量標準：Bp，B，BX，BM[12].其中，Bp是所能獲得的最小點誤差，B是權(quán)重向量取值范圍為RN時的誤差，BX是權(quán)重向量取值為凸空間時的誤差，BM是表現(xiàn)最好的模型M的誤差.以上4種度量標準之間在樣本不是太少的情況下，恒存在關(guān)系Bp≤B≤BX≤BM.

3.2 算法性能評測與分析

進行算法性能評測時，因為RR算法既可以應(yīng)用在多個站點，也可以應(yīng)用在單個站點.所以，分別對多站點的集合預測和單站點的集合預測進行了比對分析.最后，針對OCF算法、RR算法的預報性能，選取了幾個站點進行時間序列的分析比對.

3.2.1 程序運行環(huán)境與整體架構(gòu)

程序主要分為數(shù)據(jù)清洗模塊和算法應(yīng)用兩大部分，所用的編程工具是Python 2.7和MATLAB 2016b，其中數(shù)據(jù)處理時還需要以下Python庫pandas-0.19.2、numpy-1.12.0+mkl、nltk-3.2.2、scipy-0.18.1、XlsxWriter-0.9.5.Python的pandas庫中包含的數(shù)據(jù)框格式(DataFrame)非常適于進行二維的數(shù)據(jù)處理.程序數(shù)據(jù)的處理流程如圖1所示.

圖1 程序數(shù)據(jù)處理流程圖

3.2.2 多站點集合評測分析

多站點集合評測分析即是將多個站點的數(shù)據(jù)放到同一個數(shù)據(jù)集中，然后對這個數(shù)據(jù)集進行計算和處理.理想情況下，多個站點的選取最好是在同一區(qū)域中，這樣可以對某個地區(qū)進行集合預測，理論上如果監(jiān)測站點比較多，其他相鄰地區(qū)可以通過插值的方法得到預測值，對區(qū)域的預測會更有針對性.

數(shù)據(jù)清洗采用兩種力度：一種不除去重復模式值，一種除去重復模式值.以s代表參與計算的站點數(shù)，T代表日期數(shù)，BOCF和BRR分別代表OCF算法和RR算法對該數(shù)據(jù)集給出預測的預測矩陣的RMSE值.實驗中OCF算法的窗口期為7 d，計算OCF算法和RR算法RMSE值的評測時段均為第8天至第T天(第1天權(quán)重都為0，無法給出預測值).對此，進行了如下幾組測試.

(1) 數(shù)據(jù)初始規(guī)模s=30，T=90；不去重清洗后規(guī)模為s=11，T=89.此時得到表1中所列數(shù)據(jù).

去重后數(shù)據(jù)規(guī)模為s=7，T=89.此時得到表2中所列數(shù)據(jù).

(2) 數(shù)據(jù)初始規(guī)模s=90，T=90；不去重清洗后規(guī)模為s=40，T=89.此時得到表3中所列數(shù)據(jù).

表3 s=90，T=90數(shù)據(jù)RMSE值計算(不去重)

去重后數(shù)據(jù)規(guī)模為s=22，T=89.此時得到表4中所列數(shù)據(jù).

表4 s=90，T=90數(shù)據(jù)RMSE值計算(去重)

(3) 數(shù)據(jù)初始規(guī)模s=120，T=120；不去重清洗后規(guī)模為s=50，T=119.此時得到表5中所列數(shù)據(jù).

表5 s=120，T=120數(shù)據(jù)RMSE值計算(不去重)

去重后數(shù)據(jù)規(guī)模為s=22，T=119.此時得到表6中所列數(shù)據(jù).

表6 s=120，T=120數(shù)據(jù)RMSE值計算(去重)

通過以上實驗結(jié)果分析可以得到：對于較長時間維度(90 d以上)，無論數(shù)據(jù)集是否去重，RR算法性能比OCF算法性能都要好很多，相比最優(yōu)模式值BM,其優(yōu)勢也很明顯，甚至比最優(yōu)常數(shù)線性組合B的性能還要好.實際上，如果在較短時間維度上(如30 d)，RR算法的表現(xiàn)或許會比OCF算法差.根據(jù)集成學習的理論，這種情況往往是由訓練數(shù)據(jù)過少所導致.

3.2.3 單站點效果評測分析

RR算法和OCF算法都可以應(yīng)用于單站點.對于某些站點，將其作為一個集合得到的結(jié)果更優(yōu)還是單個站點獨立應(yīng)用該模式得到的方法更優(yōu)呢？為了探討這個問題，進行如下測試實驗.

隨機選取10個站點的2015年9月1日起連續(xù)120 d的數(shù)據(jù)，進行清洗后，數(shù)據(jù)規(guī)模為s=10，T=119.評測RMSE值的時段均為第16天至第T天(對于RR算法，第1天權(quán)重都為0，無法給出預測值；而OCF算法需要相應(yīng)的窗口期).

(1) 將多個站點作為1個集合應(yīng)用RR算法和OCF算法(窗口期W分別取7和15)，可以得到表7中所列數(shù)據(jù).

表7 多站點集合應(yīng)用RR算法和OCF算法的RMSE值

(2) 對這些站點分別應(yīng)用RR算法和OCF算法，則可以得到表8中所列數(shù)據(jù).

表8 單站點分別應(yīng)用RR算法和OCF算法的RMSE值

經(jīng)過分析可以發(fā)現(xiàn)：

(1) OCF算法(窗口期為7 d)的性能比OCF算法(窗口期為15 d)的好.

(2) 從整體上看，RR算法性能比OCF算法要好.多站點集合應(yīng)用時RR算法比OCF算法性能優(yōu)勢明顯，單站點應(yīng)用時，除第9號單站外，其他單站RR算法的應(yīng)用效果相比OCF算法均有較大優(yōu)勢.

(3) RR算法可以在單站點上使用，多站點集合應(yīng)用時計算出的RMSE值處于最優(yōu)單站和最差單站之間(最差單站RMSE值比集合應(yīng)用時要差).就所選樣本而言，大多數(shù)站點單獨應(yīng)用RR算法時，性能都較好.

(4) 對于RR算法中參數(shù)λ的取值，其必然遠離0，然而對于不同的站點、輸入數(shù)據(jù)或不同的評測時長，λ值通常不相同.所以，應(yīng)該根據(jù)實際情況確定該值的大小.

3.2.4 時間序列比對分析

為了對OCF算法和RR算法的性能有一個更加直觀的表示，選取了幾個站點，對模式數(shù)據(jù)和觀測數(shù)據(jù)進行時間序列的效果比對和分析，結(jié)果如下.

(1) 隨機抽取站點A，從2015年9月1日起，90 d時間維度內(nèi)(實際清洗后為89 d)針對OCF算法、RR算法、實際觀測值進行時間序列比對分析，結(jié)果如圖2所示.

圖2 站點A 90 d預測值與觀測值對照圖

(2) 隨機抽取站點B，從2015年9月1日起，90 d時間維度內(nèi)(實際清洗后為89 d)針對OCF算法、RR算法、實際觀測值進行時間序列比對分析，結(jié)果如圖3所示.

圖3 站點B 90 d預測值與觀測值對照圖

(3) 隨機抽取站點C，從2015年9月1日起，120 d時間維度內(nèi)(實際清洗后為119 d)針對OCF算法、RR算法、實際觀測值進行時間序列比對分析，結(jié)果如圖4所示.

圖4 站點C 120 d預測值與觀測值對照圖

(4) 隨機抽取站點D，從2015年9月1日起，120 d時間維度內(nèi)(實際清洗后為119 d)針對OCF算法、RR算法、實際觀測值進行時間序列比對分析，結(jié)果如圖5所示.

圖5 站點D 120 d預測值與觀測值對照圖

(5) 隨機抽取站點E，在2016年1月1日至2016年12月31日(實際清洗后為351 d)針對OCF算法、RR算法、實際觀測值進行時間序列比對分析，結(jié)果如圖6所示.

圖6 站點E 2016年預測值與觀測值對照圖

(6) 隨機抽取站點E，對2016年1月1日至2016年12月31日(實際清洗后為351 d)針對RR算法、4個模式值、觀測值進行時間序列比對分析，結(jié)果如圖7所示.

圖7 站點E 2016年RR算法值、模式值與觀測值對照圖

通過以上時間序列的比對，可以看出，RR算法比OCF算法更加貼近于真實值，對于臭氧濃度趨勢的變化適應(yīng)性更強.同時，也印證了RR算法的RMSE值為何相比OCF小很多.然而，RR算法預測的精確度依賴于4個子模式，即對于某種極端值，若子模式都沒有給出合適范圍內(nèi)的取值，則RR算法也很難直接給出這種極端值的合理預測.

3.2 算法的性能和優(yōu)勢討論

在氣象上，通常認為某個時次(某天)的誤差是由一個較為穩(wěn)定的系統(tǒng)偏差和一個確定性的擾動誤差組成.OCF算法的偏差校正其實是對系統(tǒng)偏差的校正，系統(tǒng)偏差可以通過計算一段時間的誤差平均獲得，為了得到一個魯棒性較好的平均誤差(不易受到極端值的影響)，OCF算法往往選用分位數(shù)來獲得系統(tǒng)偏差.但是，實際應(yīng)用中，OCF算法并不能很好地滿足假設(shè)條件，因此，只能從某種程度上減少極值帶來的影響.

RR算法在算法性能上，相比OCF算法和傳統(tǒng)模式(在這里作為集成模式的子模式)都有很大的提升(見圖8所示)，而且其算法僅需要一個參數(shù)，避免了很多復雜的參數(shù)條件問題.盡管該參數(shù)在指定過程中可能存在一定的不確定性，但是通過遍歷等方式并不難找到.而且RR算法的性能可以有理論上的保證，因此可以應(yīng)用的范圍也比較廣.

RR算法通?？梢园哑渌A測模式集成起來，提高整體性能.其優(yōu)勢可以歸結(jié)為：對于回歸問題，不需要構(gòu)造擬合精度和預測能力都很好的回歸算法，只要基預測器比隨機猜測略好即可；它具有很好的通用性和魯棒性，通?？梢詰?yīng)用于任何基礎(chǔ)回歸算法，神經(jīng)網(wǎng)絡(luò)、決策樹、支持向量機、線性回歸、數(shù)值模式等方法都可以作為基學習器進行集成，而且理論上可以保證集成效果比單個學習器更好，不容易出現(xiàn)過擬合問題.當然，參與集成的個體數(shù)并不是越多越好，更多的個體需要更大的計算和存儲開銷，個體間的差異更難以獲得.

圖8 數(shù)據(jù)集RMSE對比圖(不去重清洗/去重清洗)

4 結(jié)束語

論文實現(xiàn)了基于Boosting的集合預測方法嶺回歸Ridge Regression (RR)算法，針對算法中參數(shù)的特點進行了分析，并給出了調(diào)整參數(shù)的方案.進行了算法性能的測試和評價.根據(jù)我國站點觀測值和模式預報值，對OCF算法和RR算法分別給出相應(yīng)的集合預報值，并通過均方根誤差(RMSE)和時間序列進行性能評估.發(fā)現(xiàn)通常情況下，與已有的預報模式相比，RR算法具有更高的準確度和更好的穩(wěn)定性.

基于機器學習技術(shù)的Ridge Regression算法是一個元模式級別的算法，它可以將多種模式集成在一起，起到“博采眾長”的作用，它的集成性能和目標預測的空氣污染物沒有相關(guān)性，即也可以進行其他氣象指標的預測應(yīng)用.同時，它的預測準確性與這些基學習器緊密相關(guān).現(xiàn)階段，傳統(tǒng)氣象預報模式大都基于數(shù)值模式和物理過程的模擬，如果能夠通過深度學習等機器學習技術(shù)，探討出一個更強的分類器，無疑將會對集合預測的精度帶來很大提升.

進一步的研究如果能夠從算法層面更好地解決權(quán)重更新和參數(shù)選取問題，將對更大規(guī)模的數(shù)據(jù)設(shè)計高效的并行算法，為臭氧等大氣污染物的預測帶來新的機遇與挑戰(zhàn).當預測精度達到一定程度后，可以用來反推污染源位置等信息，跟蹤生產(chǎn)方式帶來的影響等，這將為各種空氣污染物的治理提供解決思路.