謝 磊,鐵治欣,宋飛揚,丁成富

1(浙江理工大學(xué) 信息學(xué)院,杭州 310018)

2(聚光科技 (杭州)股份有限公司,杭州 310052)

SO2作為大氣中主要污染物之一,既能導(dǎo)致空氣質(zhì)量變差,又能在空氣中迅速氧化生成SO3進而通過異質(zhì)核化過程形成液態(tài)氣溶膠[1],降低了大氣的能見度.同時,高濃度的SO2也是形成酸雨的罪魁禍?zhǔn)譡2],嚴重影響水環(huán)境以及農(nóng)作物的產(chǎn)量,給人們的身心以及日常活動帶來嚴重的影響.所以,實現(xiàn)對大氣中SO2的精準(zhǔn)預(yù)測,從而達到對空氣質(zhì)量的精準(zhǔn)預(yù)警預(yù)報變得尤為重要.目前,在空氣質(zhì)量預(yù)報領(lǐng)域展現(xiàn)出良好應(yīng)用前景的預(yù)報模式包括:CMAQ、CAMx、WRFCHEM、NAQPMS等[3].但是由于各模式本身的不確定因素較大,在一定程度上影響了空氣質(zhì)量的預(yù)報結(jié)果,導(dǎo)致對SO2的預(yù)測不精準(zhǔn),為了提高對大氣中各污染物濃度的預(yù)測準(zhǔn)確度,以集合的形式把多個空氣質(zhì)量預(yù)報模式的預(yù)報結(jié)果集合在一起的預(yù)報方法已經(jīng)成為提高預(yù)測大氣污染物濃度的有效決策和方法[4].

目前,運用在大氣污染物集合預(yù)報中較為廣泛的方法包括以下幾種:a)算術(shù)平均法.它具有算法簡單且能夠很好地反映出預(yù)測結(jié)果的整體趨勢等優(yōu)點.如:王自發(fā)等[5]基于多個空氣質(zhì)量預(yù)測模式(CMAQ、CAMx、NAQPMS)構(gòu)建出業(yè)務(wù)化的北京空氣質(zhì)量集合預(yù)報系統(tǒng)(EMS-Beijing),該系統(tǒng)預(yù)測結(jié)果主要依賴于各模式的算術(shù)平均值.但是僅僅利用算術(shù)平均法進行集合預(yù)報獲得的結(jié)果值極易受到極個別異常值的影響.Mallet等[6]在不同的初始條件下,互相組合構(gòu)建了幾十種空氣質(zhì)量預(yù)報模式系統(tǒng),發(fā)現(xiàn)只是簡單的平均集合方法并不能顯著地改進預(yù)報結(jié)果,并且平均集合方法的使用范圍總存在著一定的局限性,需要引入一些其它的方法才能改進預(yù)報性能.b)多元線性回歸法.該方法通過研究多個自變量和因變量的線性關(guān)系,計算出一個擬合度較高的回歸預(yù)測模型,能夠很好地度量各個因素之間的相關(guān)程度.如:黃思等[7]利用多模式集合和多元線性回歸改進北京PM10預(yù)報,通過將歷史觀測信息納入到北京空氣質(zhì)量集合預(yù)報系統(tǒng)(EMSBeijing)預(yù)報體系,并對EMS-Beijing三個集合成員的預(yù)報結(jié)果進行集成,能夠很好的分析出影響集成預(yù)報的影響因子,提高了預(yù)測的準(zhǔn)確率.c)動態(tài)權(quán)重分配法.這種方法通過分析在一段時間內(nèi)各模式預(yù)測值的偏差率動態(tài)調(diào)整集合模式的權(quán)值,有效地提高了集合模式的預(yù)測準(zhǔn)確率.如:姚文強等[8]針對杭州市的空氣質(zhì)量情況,提出了一種動態(tài)權(quán)重更新模型并與已有的平均預(yù)測模型進行比較,結(jié)果表明比簡單平均集成效果更好.但是通過單個預(yù)測模型的集成方法的適用性存在一定的局限性,且預(yù)測準(zhǔn)確率不穩(wěn)定[9].張春萍等[10]在進行區(qū)域物流量的預(yù)測當(dāng)中,基于多個單項預(yù)測模型的預(yù)測結(jié)果,采用最優(yōu)定權(quán)組合模型以誤差平方和最小為原則求得各單項預(yù)測模型的最優(yōu)權(quán)值,有效地解決了單項預(yù)測模型準(zhǔn)確率不穩(wěn)定,適用性不高等問題,查閱相關(guān)資料得知有關(guān)最優(yōu)定權(quán)組合預(yù)測模型在空氣質(zhì)量預(yù)測方面鮮有介紹.

因此,本文通過采用三種國內(nèi)應(yīng)用成熟、業(yè)務(wù)化程度高的空氣質(zhì)量預(yù)測模式(CMAQ、CAMx、WRFCHEM),結(jié)合云南省的地勢、氣候、污染源等特點,提出了一種基于最優(yōu)定權(quán)組合預(yù)測模型的集合預(yù)測方法將上述三種預(yù)測模式的結(jié)果進行集成,并通過誤差分析驗證所提方法的有效性.

1 模式系統(tǒng)與數(shù)據(jù)介紹

1.1 系統(tǒng)的介紹

云南省省級空氣質(zhì)量預(yù)報預(yù)警平臺所采用的多模式集合預(yù)報系統(tǒng)成員分別是由美國國家海洋和大氣管理局(NOAA)預(yù)報系統(tǒng)實驗室開發(fā)的WRF-CHEM[11],美國國家環(huán)境保護局研發(fā)的CMAQ[12]系統(tǒng)以及美國環(huán)境技術(shù)公司Environ維護并研發(fā)的CAMx[13]組合而成.上述三種模型被廣泛的應(yīng)用于國內(nèi)外的空氣質(zhì)量模擬領(lǐng)域,但各個模式輸出的格式各有特點,在定制化的過程中,因為變量名、變量組別的不同,為集合預(yù)報帶來了極大的挑戰(zhàn).為了實現(xiàn)模式的統(tǒng)一化管理與展示,參考國家環(huán)保部發(fā)布的《環(huán)境空氣質(zhì)量指數(shù)(AQI)技術(shù)規(guī)定》,將三種不同的模式集成到同一個系統(tǒng)平臺上,氣象場采用中尺度WRF模式,初始場采用GFS數(shù)據(jù)集.云南省多模式集合預(yù)報系統(tǒng)框架如圖1所示.

首先基于不同污染源的排放特征,通過SMOKE源排放模型對排放源清單進行處理,為各單項空氣質(zhì)量模式提供所需要的大氣排放源清單的時空分布數(shù)據(jù).然后利用氣象模式WRF對GFS數(shù)據(jù)集進行分析計算,并按照各單項空氣質(zhì)量模式的運行需求將上述計算結(jié)果轉(zhuǎn)換成其所需要的氣象文件.最后將氣象模式和源排放模型的計算結(jié)果作為各單項空氣質(zhì)量模式的輸入計算出各模式對應(yīng)的SO2預(yù)測值,并分別采用最優(yōu)定權(quán)組合預(yù)測法、多元線性回歸法和動態(tài)權(quán)重更新法對上述各單項空氣質(zhì)量模式的SO2預(yù)測值進行集成,從而實現(xiàn)多模式集合預(yù)報.這一系統(tǒng)已經(jīng)應(yīng)用于云南省環(huán)境檢測中心,為云南省的空氣質(zhì)量預(yù)報提供更加準(zhǔn)確的預(yù)報數(shù)據(jù).

圖1 云南省多模式集合預(yù)報系統(tǒng)框架

1.2 網(wǎng)格嵌套設(shè)置

本文重點針對云南省環(huán)境監(jiān)測中心業(yè)務(wù)化運行的云南省省級環(huán)境空氣質(zhì)量預(yù)報預(yù)警業(yè)務(wù)平臺進行分析與集成研究.為了便于系統(tǒng)采用的三種不同模式的集成以及減小因為不同的區(qū)域設(shè)置而帶來的誤差,本系統(tǒng)采用統(tǒng)一的區(qū)域設(shè)置、統(tǒng)一的排放清單、統(tǒng)一的氣象場驅(qū)動,以云南省為中心的三重網(wǎng)格設(shè)置.第一重網(wǎng)格包括中國中西部地區(qū),范圍如圖2中的D1所示,網(wǎng)格精度為 27 km,網(wǎng)格數(shù) (東西×南北)為 141×134,第二重網(wǎng)格包括中國西南地區(qū),范圍如圖2中的D2所示,網(wǎng)格精度 9 km,網(wǎng)格數(shù)為 219×219,第三重網(wǎng)格包括云南全省,范圍如圖2中的D3所示,網(wǎng)格精度為3 km,網(wǎng)格數(shù)為327×327.

圖2 云南省多模式集合預(yù)報系統(tǒng)網(wǎng)格設(shè)置簡圖

1.3 運行環(huán)境

云南省環(huán)境監(jiān)測中心預(yù)報預(yù)警業(yè)務(wù)平臺采用多模式集合預(yù)報系統(tǒng)的環(huán)境為 Red Hat Linux 7.3,編譯環(huán)境為 Intel Fortran,并行環(huán)境為 MPI和 OpenMP 混合并行,運行所用核數(shù)為48核/96線程,同時安裝了Intel C/C++編譯器、GNU C/C++編譯器.

1.4 實驗數(shù)據(jù)

本文所采用的數(shù)據(jù)均來自于云南省省級環(huán)境空氣質(zhì)量預(yù)報預(yù)警業(yè)務(wù)平臺數(shù)據(jù)庫,為保證數(shù)據(jù)的有效性,該數(shù)據(jù)庫數(shù)據(jù)均經(jīng)過自動審核系統(tǒng)的過濾和篩選.其中觀測數(shù)據(jù)分別來源于云南省蒙自、楚雄、昭通三個區(qū)域站點在2018年1月至6月期間的SO2日均實測值.模式輸出數(shù)據(jù)來自模式第三重網(wǎng)格的預(yù)報結(jié)果,根據(jù)我國執(zhí)行的空氣污染指數(shù)標(biāo)準(zhǔn),利用當(dāng)天的12時至第二天的11時的SO2預(yù)報值計算日均值.

2 集成方法和評價指標(biāo)

2.1 本文的最優(yōu)定權(quán)組合預(yù)測模型

組合預(yù)測最早是由Bates和Granger提出來的[14],其基本思想是即使某一單項預(yù)測方法的預(yù)測精度很差,我們也不能隨意丟棄任何一種預(yù)測結(jié)果,因為這種做法極有可能造成某些有用信息的丟失,把它和某一種或幾種預(yù)測精度較高的方法進行組合,完全有可能提高系統(tǒng)的整體預(yù)測能力,比較科學(xué)的做法是通過某種方法將各單項預(yù)測方法進行組合集成,綜合利用各單項預(yù)測方法的優(yōu)勢,從而獲得優(yōu)于各單項預(yù)測方法的預(yù)測能力.

組合預(yù)測的關(guān)鍵[15]是如何確定各組合成員的權(quán)重系數(shù),對于多模式空氣質(zhì)量數(shù)值模式集合來說,傳統(tǒng)的方法是主觀賦權(quán)法,通過相關(guān)領(lǐng)域的專家根據(jù)經(jīng)驗進行主觀判斷得到相應(yīng)的權(quán)重系數(shù),這種方法的主觀因素較多,隨意性較大,不同的專家給出的結(jié)果不同.還有一種是普通的平均定權(quán)法,雖然方法簡單,但是如果某一組合預(yù)測成員的誤差極大,通過平均定權(quán)后很容易影響整個系統(tǒng)的預(yù)測精度.而最優(yōu)定權(quán)組合法基于誤差平方和最小原則求得最優(yōu)加權(quán)系數(shù),它綜合了各種單項預(yù)測模型的優(yōu)勢,充分利用各種單項模型提供的有用信息達到提高預(yù)測準(zhǔn)確度的效果,不僅能克服上述缺點,而且還能在一定程度上分散組合預(yù)測模型中選擇不當(dāng)?shù)娘L(fēng)險.

利用最優(yōu)定權(quán)組合法集成的思路如下:通過應(yīng)用多個單項空氣質(zhì)量數(shù)值模式(WRF-CHEM、CMAQ、CAMx)對SO2的濃度進行預(yù)測,并計算出相應(yīng)的單項模式預(yù)測誤差,然后將各單項預(yù)測誤差以不同的權(quán)值建立出組合預(yù)測誤差方程,并以誤差平方和最小為原則建立目標(biāo)函數(shù)對組合預(yù)測誤差方程進行優(yōu)化,最后在一定的約束條件下求出目標(biāo)函數(shù)的最優(yōu)解.

具體建模步驟如下:假設(shè)SO2的觀測值序列為Y={yt,t= 1,2,… ,n} ,其中yt為第t時刻的實際觀測數(shù)據(jù).利用m種單項預(yù)測模型對SO2濃度進行預(yù)測,設(shè)第i(i=1,2,… ,m)個單項預(yù)測模型在第t時刻的預(yù)測值為,且第i個單項預(yù)測模型的權(quán)值為ki,并且滿足,則第t時刻組合預(yù)測值可表示為:

令第i個單項模型在t時刻的單項預(yù)測誤差為eit,則eit可表示為:

令et為組合預(yù)測模型在第t時刻的組合預(yù)測誤差,則et可表示為:

令J為組合預(yù)測模型的誤差平方和,則有:

則式(4)的矩陣表達式如下所示:

令R=(1,1,···,1)T,則可以表示為:

由此可知我們要做的就是在式(6)的約束條件下,求出最優(yōu)加權(quán)向量K,使得此時的誤差平方和J最小,即:

利用拉格朗日乘子λ對式(7)進行求解,則J可以表示為[16]:

若想求得J的極值,則J對K的一階偏導(dǎo)必為零[16],則:

將式(9)兩邊左乘E-1可得:

將式(10)兩邊左乘RT可得:

將式(6)帶入式(11)即可得到 λ的值:

將式(12)帶入式(10)即可得到K的值:

又由式(6)可得:

從而得出結(jié)論,由式(13)求得的最優(yōu)權(quán)重滿足式(6)的約束條件,即表明式(13)求得的K就是最優(yōu)加權(quán)系數(shù)向量,且此時的J值就是極小值.

最后將求出的最優(yōu)權(quán)重K代入式(1),即可求出t時刻的最優(yōu)預(yù)測值.

2.2 多元線性回歸法

多模式集合預(yù)報系統(tǒng)的預(yù)報結(jié)果往往受到多個單項預(yù)報成員的影響,由于各單項預(yù)報模式的預(yù)報性能不同,不同模式對SO2的預(yù)測準(zhǔn)確程度不同,所以為了分析各單項預(yù)報模式之間的相關(guān)性,采用多元線性回歸的方法[7]對各模式的預(yù)報結(jié)果進行集成預(yù)報,通過大量的樣本值做回歸分析,建立擬合度較高的回歸方程.

2.3 動態(tài)權(quán)重更新法

動態(tài)權(quán)重更新[8]通過評估一段時間內(nèi)單個預(yù)報模式的預(yù)測值和實際觀測值的偏差情況,從而可以動態(tài)調(diào)整各單項預(yù)報模式在集成預(yù)報中的權(quán)重,其優(yōu)點在于各單項預(yù)報模式的權(quán)重分配依賴于實際的監(jiān)測數(shù)據(jù),具體操作步驟如下:

設(shè)第i個數(shù)值預(yù)報模式第j天的SO2預(yù)測值為,第j天的實測日數(shù)據(jù)為Yj,則該模式的偏差率為:

然后根據(jù)第i個模式每天的偏差率計算出總共m天內(nèi)的平均偏差率Rim.定義一個過渡因子Vi,且該過渡因子與權(quán)重系數(shù)呈線性關(guān)系,與預(yù)測偏差率呈反比,則Vi可表示為:

設(shè)單個模式的權(quán)重為wi(i= 1,2,…,n),且滿足,則wi可表示為:

最后將求得的權(quán)重和相對應(yīng)的模式預(yù)報結(jié)果進行加權(quán)計算即可得到多模式集成預(yù)報結(jié)果.

2.4 評估方法

根據(jù)統(tǒng)計學(xué)的相關(guān)性原理,本文通過使用誤差平方和 (Sum of Squares for Error,SSE)以及均方百分比誤差 (Mean Square Percentage Error,MSPE)兩個評價指標(biāo)對各預(yù)測模型的準(zhǔn)確度進行評估分析.

(1)誤差平方和:

(2)均方百分比誤差:

其中,n表示實驗樣本總個數(shù),為第t時刻單項預(yù)測方法的預(yù)測值,xt為第t時刻的實測值.

SSE表征的是由于模式的輸入、排放清單采集不合理以及模式自身的一些物理和化學(xué)參數(shù)導(dǎo)致預(yù)報結(jié)果與實際觀測值的偏差.其值越大,表明與測定值之間的差異越大.MSPE則表征各模式的平均誤差水平.其值越小,表明模式的預(yù)測值與實際觀測值的誤差值越小,即預(yù)測的偏差值更小.

3 實驗結(jié)果

3.1 單項模式預(yù)測評估

為了了解單個空氣質(zhì)量模式的特點以及它在不同站點對SO2的預(yù)測能力,本文收集了2018年6月份云南省蒙自、楚雄、昭通三個站點各單項空氣質(zhì)量預(yù)報模式(CMAQ、CAMx、WRF-CHEM)的SO2預(yù)測值和日均實測值.通過對比不同區(qū)域的各單項預(yù)報模式預(yù)測結(jié)果,利用SSE與MSPE兩項指標(biāo)評估多模式集合系統(tǒng)中三個空氣質(zhì)量模式對不同站點SO2的預(yù)測能力,結(jié)果如表1所示.

觀察表1中數(shù)據(jù)發(fā)現(xiàn),在蒙自站點的誤差分析中,CMAQ模式的SSE以及MSPE在三種模式誤差對比中最小,說明該模式的預(yù)測能力最好,CAMx模式的預(yù)測能力次之,WRF-CHEM模式的預(yù)測能力最差.而在楚雄站點的誤差分析中,對SO2預(yù)測較為精準(zhǔn)的模式是WRF-CHEM和CAMx,CMAQ模式的預(yù)測能力最差.在昭通站點的誤差分析中,WRF-CHEM和CMAQ模式的預(yù)測能力較為精準(zhǔn),CAMx的預(yù)測能力最差.

由此可知,對于不同的區(qū)域,不同的模式預(yù)測能力相差較大,每種模式在不同區(qū)域的優(yōu)點各有不同.通過進一步對比不同區(qū)域的三種模式的預(yù)測結(jié)果可得,沒有任何一種模式在不同區(qū)域的預(yù)測結(jié)果完全優(yōu)于其它兩個模式,從而在一定程度上證明了單項空氣質(zhì)量模式對大氣污染物濃度的預(yù)測仍存在一定的局限性.

表1 各單項模式在不同區(qū)域的誤差分析

3.2 最優(yōu)定權(quán)組合法集成預(yù)測評估

為了充分利用各單項空氣質(zhì)量預(yù)測模式的優(yōu)勢,驗證本文所提出的最優(yōu)組合預(yù)測模型(OWCF)的有效性,選取多元線性回歸法(MLR)和動態(tài)權(quán)重更新法(DWA)與本文提出的最優(yōu)組合預(yù)測法分別對各單項空氣質(zhì)量預(yù)報模式(CMAQ、CAMx、WRF-CHEM)的預(yù)報結(jié)果進行集成對比實驗.即在相同的實驗條件下,分別利用蒙自、楚雄、昭通三個站點2018年1至5月份的SO2實測值和上述模式預(yù)報數(shù)據(jù)構(gòu)建出預(yù)測模型,對各站點2018年6月份的SO2濃度值進行預(yù)測,并與各站點六月份的實測數(shù)據(jù)(REAL)進行對比驗證,最后通過繪制對比圖的方式將各方法的預(yù)測結(jié)果與實測結(jié)果進行對比分析,實驗結(jié)果如圖3-圖11所示.

圖3 MLR 預(yù)測結(jié)果圖(蒙自)

圖4 DWA 預(yù)測結(jié)果圖(蒙自)

圖5 OWCF 預(yù)測結(jié)果圖(蒙自)

圖6 MLR 預(yù)測結(jié)果圖(楚雄)

由圖3-圖11中不同站點的預(yù)測結(jié)果對比曲線可以看出,在大多數(shù)時間內(nèi),三種方法的預(yù)測結(jié)果均與實測結(jié)果(REAL)的擬合曲線趨勢一致.但是相較于本文所提的OWCF,MLR和DWA所預(yù)測的結(jié)果準(zhǔn)確率不高,偏差較大,而OWCF的結(jié)果值更加接近于實測值,預(yù)測準(zhǔn)確率更高.

圖7 DWA 預(yù)測結(jié)果圖(楚雄)

圖8 OWCF 預(yù)測結(jié)果圖(楚雄)

圖9 MLR 預(yù)測結(jié)果圖(昭通)

為了定量地分析最優(yōu)定權(quán)組合預(yù)測模型的特點,分別運用SSE與MSPE對三種預(yù)測方法的預(yù)測精度實行評估檢驗,其結(jié)果如表2所示.

表2表示通過采用式(20)和式(21)對三種預(yù)測方法在不同區(qū)域的預(yù)測結(jié)果進行評估分析,由表中數(shù)據(jù)可知,在蒙自站點中,MLR和DWA的SSE分別為179.82、218.05,本文所提的OWCF的SSE與以上兩種方法的SSE相比分別降低了143.65、181.88.從MSPE的角度來看,MLR、DWA的均方百分比誤差分別為 0.047、0.054,OWCF 的MSPE為 0.022,與以上兩種方法相比分別降低了0.025、0.032個百分點.同樣,在楚雄站點和昭通站點的數(shù)據(jù)中,本文所提的OWCF預(yù)測法的預(yù)測精度都存在不同程度的提升,SSE、MSPE均最小.

圖10 DWA 預(yù)測結(jié)果圖(昭通)

圖11 OWCF 預(yù)測結(jié)果圖(昭通)

表2 三種預(yù)測方法在不同區(qū)域的誤差分析

且通過對比表1中單項模式的誤差可知,在蒙自站點中,預(yù)測能力最好的單項模式是CMAQ,其中SSE為254,MSPE為0.056.而通過采用本文所提的OWCF法將各模式結(jié)果進行集成所獲得的SSE為36.17,MSPE為0.022.說明本文所提方法的預(yù)測精度即使和單項預(yù)測模式中預(yù)測精度最優(yōu)的模式相比任然存在一定的優(yōu)勢.同樣,在楚雄和昭通站點中,OWCF和各站點單項最優(yōu)的預(yù)測模式相比,OWCF的SSE和MSPE均最低.

由以上分析可知,通過與單項預(yù)測精度最優(yōu)的模式以及采用不同方式的集成方法預(yù)測效果相比較,本文所提的OWCF有效的提升了SO2的預(yù)測準(zhǔn)確度,從而證明了本文所提方法的可行性.

4 結(jié)語

本文基于多個空氣質(zhì)量預(yù)測模式(WRF-CHEM、CMAQ、CAMx),以誤差平方和最小為原則,提出了最優(yōu)定權(quán)組合預(yù)測模型.采用2018年1至5月份云南省三個監(jiān)測站點(蒙自、楚雄、昭通)的SO2濃度實測值和多個單項空氣質(zhì)量模式的預(yù)測值作為實驗數(shù)據(jù),與現(xiàn)有的多元線性回歸法、動態(tài)權(quán)重分配法等預(yù)測方法對前述三個站點2018年6月份的SO2濃度進行預(yù)測,并對各模型的預(yù)測誤差進行評估分析.實驗結(jié)果表明最優(yōu)定權(quán)組合預(yù)測模型的預(yù)測準(zhǔn)確率最高,效果最佳,為云南省空氣質(zhì)量預(yù)警預(yù)報平臺提供一種高效的預(yù)測方法.