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

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

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

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

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

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

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

1.1 系統(tǒng)的介紹

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

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

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

1.2 網格嵌套設置

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

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

1.3 運行環(huán)境

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

1.4 實驗數據

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

2 集成方法和評價指標

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

又由式(6)可得:

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

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

2.2 多元線性回歸法

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

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

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

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

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

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

最后將求得的權重和相對應的模式預報結果進行加權計算即可得到多模式集成預報結果.

2.4 評估方法

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

(1)誤差平方和:

(2)均方百分比誤差:

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

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

3 實驗結果

3.1 單項模式預測評估

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

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

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

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

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

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

圖3 MLR 預測結果圖(蒙自)

圖4 DWA 預測結果圖(蒙自)

圖5 OWCF 預測結果圖(蒙自)

圖6 MLR 預測結果圖(楚雄)

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

圖7 DWA 預測結果圖(楚雄)

圖8 OWCF 預測結果圖(楚雄)

圖9 MLR 預測結果圖(昭通)

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

表2表示通過采用式(20)和式(21)對三種預測方法在不同區(qū)域的預測結果進行評估分析,由表中數據可知,在蒙自站點中,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個百分點.同樣,在楚雄站點和昭通站點的數據中,本文所提的OWCF預測法的預測精度都存在不同程度的提升,SSE、MSPE均最小.

圖10 DWA 預測結果圖(昭通)

圖11 OWCF 預測結果圖(昭通)

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

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

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

4 結語

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