天目山地區重污染天氣預測模型及其環流特征研究

李 云， 張 瑩， 許金萍， 蔣曉梅， 楊 聃， 梁明珠

(1.浙江省安吉縣氣象局,浙江 安吉 313300； 2.浙江省湖州市氣象局,浙江 湖州 313000)

引 言

在全國污染物減排的大背景下，氣象條件顯然已成為大氣污染濃度變化的最重要因子[1-7]。包括地面氣象觀測要素及天氣形勢場及其相關物理量在內的氣象條件，是研究地區大氣重污染事件的主要對象。然而污染物濃度變化的主要驅動因子不同研究區域之間有較大的差別。Zhang等[8]研究發現，環流形勢是北京地區污染物濃度逐日變化的主要驅動因子。楊茜等[9]探究了降水對重慶市大氣污染物濃度的影響，發現PM2.5、PM10濃度隨降水量增加逐漸降低，降低趨勢線較為明顯。馬格等[10]對鄭州市氣象因子對大氣顆粒物濃度的影響研究發現，降水、風向、濕度與大氣顆粒物濃度相關。劉雯等[11]對武漢市觀象臺2013—2016年PM2.5質量濃度變化及其與氣象因子的相關分析發現，PM2.5質量濃度變化與降水量、風速和氣溫等氣象因子明顯相關。王景云等[12]分析了2012—2015年北京市空氣質量指數變化及其與氣象要素的關系，發現整體上空氣質量指數與風速、日照時數、降水量、平均氣溫和最高氣溫呈負相關，與濕度呈正相關。李德平[13]、孫燕[14]等通過分析發現，風是重污染的主要因子之一。劉郁玨等[15]對北京市房山區大氣污染物時空分布特征及氣象影響因素分析發現，在不同季節條件下，局地氣象要素與污染天氣發生概率之間有著很好的相關關系；Pearce等[16-17]研究也指出，局地氣象條件而非天氣形勢是污染物日變化的主要驅動因子。大氣污染物濃度的變化與多個非獨立因子的氣象條件有著非線性關系。何建軍等[18]利用ANN技術，構建了氣象條件、污染物排放變化和污染物濃度的預測模型，其對NO2模擬準確度最好，其次是SO2的模擬，對PM10的模擬結果最差。宋丹等[19]通過多元線性逐步回歸和BP神經網絡方法，分季節建立空氣質量指數預報模型，夏季指數準確率近99%，冬季超過或接近80%。

為構建一個較高準確度的重污染天氣時PM10及PM2.5最優預測模型及其氣象條件特征，本文以浙江天目山為例開展研究。

天目山地區位于浙江西北部，西接皖南，經浙、皖邊境過杭嘉湖平原西緣，呈西南—東北走向。山脈主體在安吉、臨安境內，余脈還延伸至德清、余杭境內。由于浙江省內的重污染天氣西部比東部的多，持續時間長，北部又比南部的多，加之天目山山脈影響，致使天目山地區成了省內冬、春空氣污染的重災區之一。由于地處天目山山脈的安吉和臨安，都是國家級生態建設示范區，因此本文重點探討天目山地區重污染時期氣象條件與污染濃度預測模型及其氣象條件特征。

天目山地區空氣污染監測數據時間序列普遍偏短，重污染天氣數量少。李云等[20]研究發現，安吉地區大氣顆粒物濃度及氣象要素數據的多元線性回歸方程擬合度差，與顯著性要素偏少有密切關系。

據有關研究和大量的使用結果顯示，支持向量機(SVM)[21-24]是一種優秀的淺層學習方法，在小樣本訓練集上有著無可比擬的優勢，能夠得到比其他算法好很多的結果。LIBSVM是臺灣大學林智仁教授開發設計的一個快速有效的SVM模式識別與回歸方法。因此本文嘗試使用支持向量機(LIBSVM)方法，以期解決重污染時期氣象條件與污染濃度預測模型精確度問題。

1 研究方法

在針對2015年1月—2018年10月天目山地區出現的重污染天氣(AQI大于等于151)，分別建立PM2.5、PM10日數據最優預測模型，并對其環流特征進行分析。

1.1 基于LIBSVM最優預測模型構建

擇優選取天目山地區8個大氣成分監測站，以及安吉、臨安、德清及余杭氣象國家站小時數據，對天目山地區地面氣象觀測要素(包括雨量、氣溫、最高氣溫、最低氣溫、本站氣壓、海平面氣壓、濕靜力能量、最大風速、極大風速、露點溫度、溫度露點差、水汽壓、相對濕度、最小濕度、能見度、最小能見度、地面溫度、地面最低溫度)、NECP再分析(FNL)資料(包括30°-31°E、120°-121°N的1000、850、700、500 hPa形勢場及邊界層高度、最大緯向風、最大經向風、2 m溫度、地表溫度等相關物理場)及對應的PM2.5、PM10濃度進行數據預處理和日數據統計，利用LIBSVM方法,對PM2.5、PM10分別建立預測模型及參數(主要是懲罰參數c和核函數g)尋優。

建立基于LIBSVM方法預測模型及參數尋優，首先對PM2.5、PM10、地面氣象要素、FNL再分析資料分別歸一化處理，然后進行PCA降維，最后用交叉驗證(CV)、粒子群優化(PSO)和遺傳(GA)等算法分別進行預測模型參數對比尋優。

1.2 最優預測模型c和g參數尋優

關于SVM參數的優化選取，國際上沒有公認的統一方法，目前常用的是利用K-CV算法進行交叉驗證[25]。

K-fold Cross Validation(K-CV)是常見的CV方法之一。它將原始數據分成K組，將每個子集數據分別作一次驗證集，其余的K-1組子集數據作為訓練集，這樣會得到K個模型，將K個模型最終的驗證集的分類準確率的平均數作為此K-CV下分類器的性能指標。K-CV可以有效地避免過學習及欠學習狀態的發生，最后得到的結果比較具有說服性。

K-CV算法需要遍歷網格內所有參數點參數尋優，相對費時。PSO和GA等啟發式算法的參數尋優可更快找到最優解。GA是一種基于生物遺傳和進化機制的適合復雜系統優化的自適應概率優化技術，也是一種實用、高效、魯棒性強的優化算法。PSO是一種基于群體智能的演化算法，通過粒子在解空間追隨最優的例子進行搜索。

1.3 對應最優環流形勢場及物理量指標特征分析

利用構建PM2.5、PM10最優預測模型的地面氣象觀測要素、FNL再分析資料結果，對參與最優預測模型的環流形勢場及相關物理量進行k-means聚類分析，根據相關研究和本地經驗選取k=4、5、6進行對比擇優；根據擇優結果，將相同類型的環流形勢場進行合成計算，得到重污染天氣PM2.5、PM10對應的最優環流形勢場，并分析重污染天氣環流特征。

2 結果分析

2.1 重污染天氣下PM2.5日數據最優預測模型及參數尋優

根據與PM2.5的相關系數的絕對值與顯著性特點，將55個氣象條件按照順序進入模型參與預測，在訓練值占比69%，相關性與顯著性最好的27個要素參與時得到了重污染天氣下PM2.5日數據最優預測模型(表1)。由表1可知，最佳c為2.7114，g為3.6141，訓練值和測試值的R2分別達到了0.9992和0.7196。對比分析預測值與初始值發現(圖1)，不僅訓練值與趨勢都很好吻合，測試值也取得了很好的檢驗結果。因此本文構建的模型考慮要素合理，方法合適，取得了很好的預報效果。

表1 PM2.5最優預測模型參數尋優

圖1 基于支持向量機的訓練集(a)和測試集(b)PM2.5最優預測模型

2.2 重污染天氣下PM10日數據最優預測模型及參數尋優

與PM2.5一樣的方法與步驟，在訓練值占比69%，相關性與顯著性最好的24個要素參與時得到了重污染天氣下PM10日數據最優預測模型(表2)。由表2可看出，最佳c=30.3874，g=1.6613，訓練值和測試值的R2分別達到了0.9978和0.7792。對比分析預測值與初始值發現(圖2)，不僅訓練值與趨勢都較好吻合，雖然測試值的趨勢不如PM2.5的好，但明顯優于何建軍等[18]得到的R值0.67(R2值0.4489)。因此本文構建的模型考慮要素合理，方法合適，取得了更好的預報效果。

圖2 基于支持向量機的訓練集(a)和測試集(b)PM10最優預測模型

表2 PM10最優預測模型參數尋優

2.3 對應最優環流形勢場特征分析

聚類有效性的評價標準有兩種：一是外部標準，通過測量聚類結果和參考標準的一致性來評價聚類結果的優良；另一種是內部指標，用于評價同一聚類算法在不同聚類數條件下聚類結果的優良程度，通常用來確定數據集的最佳聚類數。

最佳聚類數判定的方法，對于內部指標，通常分為三種類型：基于數據集模糊劃分的指標，基于數據集樣本幾何結構的指標，基于數據集統計信息的指標。基于數據集樣本幾何結構的指標根據數據集本身和聚類結果的統計特征對聚類結果進行評估，并根據聚類結果的優劣選取最佳聚類數。本文主要介紹Calinski-Harabasz(CH)指標、Davies-Bouldin(DB)指標和silhouette(SI)指標。

(1)CH指標

CH指標通過類內離差矩陣描述緊密度、類間離差矩陣描述分離度，指標定義為

(1)

其中，n表示聚類的數目,k表示當前的類,trB(k)表示類間離差矩陣的跡,trW(k)表示類內離差矩陣的跡。CH越大,代表著類自身越緊密，類與類之間越分散，即更優的聚類結果。

(2)DB指標

DB指標描述樣本的類內散度與各聚類中心的間距，定義為

(2)

其中，K是聚類數目，Wi表示類Ci中的所有樣本到其聚類中心的平均距離，Wj表示類Ci中的所有樣本到類Cj中心的平均距離，Cij表示類Ci和Cj中心之間的距離。可以看出，DB越小,表示類與類之間的相似度越低，從而對應越佳的聚類結果。

(3)SI指標

對于D中的每個對象o,計算o與o所屬的簇內其他對象之間的平均距離a(o):

(3)

b(o)是o到不包含o的所有簇的最小平均距離

(4)

輪廓系數定義為

(5)

Si越接近1，則說明樣本i聚類越合理。

2.3.1 PM2.5環流形勢場分類

利用構建PM2.5最優預測模型的地面氣象觀測

要素、FNL再分析資料，對參與最優預測模型的環流形勢場及相關物理量進行k-means聚類分析。根據相關研究和本地經驗選取k=4、5、6進行對比擇優[26],其中k=5時,DB指標、SI指標最優，CH指標次優。YC最優污染天氣環流類型分5類(表3)，第四類最多，第五類次之，第二、三類最少。其中第一類5天，占17.2%；第二類、第三類各2天，各占6.9%；第四類11天，占37.9%；第五類9天，占31.0%(表略)。

表3 PM2.5和PM10不同k值聚類有效性指標比較

2.3.2 PM2.5環流特征

利用k-means聚類方法，最優歸類出5類，根據相關性較好的700 hPa合成環流形勢分別是第1類弱脊型(圖3a),第2類低槽南壓型(圖3b),第3類淺槽東移型(圖3c),第4、5類均屬于高壓脊控制(影響)型(圖3d和3e)，而第4類的脊線更密，脊區更深厚。

圖3 第1—5類700 hPa合成環流形勢(a)為第1類弱脊型，(b)為第2類低槽南壓型，(c)為第3類淺槽東移型，(d)(e)為第4、5類高壓脊控制(影響)型；圖中★指天目山地區

第2、3類天目山地區分別屬于槽前、槽后天氣，是5種類型中最少的一類，第2類容易出現降水，污染持續時間較短；第3類氣溫低,露點低，濕靜力能量小，大氣邊界層高度較低，不利于污染擴散；第1類弱脊型，濕度大，溫度露點差小，容易出現降水，不利于污染持續；第4、5類均屬于高壓脊控制(影響)型，天目山地區處在西北或偏西氣流控制，容易出現長時間持續污染，是天目山地區影響PM2.5濃度的最主要天氣類型。

2.3.3 PM10環流形勢場分類

對參與最優預測模型的環流形勢場及相關物理量進行k-means聚類分析，根據相關研究和本地經驗，選取k=4、5、6進行對比擇優。其中k=6時，DB指標、CH指標最優。因此最優污染天氣環流類型分5類(表3)，第1類最多，第2、3類次之，第6類最少。其中第1類9天，占31.0%；第2類、第3類各6天，各占20.7%；第4類5天，占17.2%；第5類、第6類分別為2天、1天，合計10.3%(表略)。

2.3.4 PM10環流特征

利用k-means聚類方法最優歸類出的6類，根據相關性較好的700 hPa環流進行分析，對應500 hPa合成環流形勢分別是第1、3、4類均屬于高壓脊控制(影響)型(圖4a、4c和4d)，而第4類的脊區更深厚、寬廣；第2類為低槽南壓型(圖4b)；第5類低渦東移型(圖4e)；第6類低渦南壓型(圖4f)。

圖4 6類700 hPa合成環流形勢(b)為第2類低槽南壓型，(e)為第5類低渦東移型，(f)為第6類低渦南壓型，(a)(c)(d)為第1、3、4類高壓脊控制(影響)型；圖中★指天目山地區

第5、6類天目山地區分別屬于槽(渦)后天氣，是6種類型中最少的兩類，容易出現降水，污染持續時間較短；第2類氣溫低，露點低，濕靜力能量較小，大氣邊界層高度較低，不利于污染擴散；第1、3、4類均屬于高壓脊控制(影響)型，天目山地區處在西北或偏西氣流控制，容易出現長時間持續污染，是天目山地區影響PM10濃度的最主要天氣類型。

3 結 論

針對2015年1月-2018年10月天目山地區出現重污染天氣(AQI大于等于151)，建立基于LIBSVM方法的預測模型及參數尋優；對參與最優預測模型的環流形勢場及相關物理量進行k-means聚類分析，得到重污染天氣PM2.5、PM10對應最優環流形勢場并分析重污染天氣環流特征，得到如下結論。

(1)根據與PM2.5的相關系數的絕對值及顯著性特點，將55個氣象條件按照順序進入模型參與預測，在訓練值占比69%，相關性與顯著性最好的27個要素參與時得到了重污染天氣下PM2.5日數據最優預測模型，訓練值和測試值的R2分別達到了0.9992和0.7196，訓練值與趨勢都能很好吻合，測試值也取得了很好的檢驗結果。因此本文構建的模型考慮要素合理，方法合適，取得了很好的預報效果。

(2)由相關性與顯著性最好的24個要素參與，得到了重污染天氣下PM10日數據最優預測模型，雖然測試值的趨勢不如PM2.5的好，訓練值和測試值的R2分別達到了0.9978和0.7792，訓練值與趨勢吻合較好。因此本文構建的PM10模型考慮要素合理，方法合適，也取得了更好的預報效果。

(3)PM2.5最優預測模型對應最優5類700 hPa合成環流形勢中第4、5類均屬于高壓脊控制(影響)型，容易出現長時間持續污染，是天目山地區影響PM2.5濃度的最主要天氣類型。PM10最優預測模型對應最優6類700 hPa合成環流形勢中第1、3、4類均屬于高壓脊控制(影響)型，容易出現長時間持續污染，是天目山地區影響PM10污染的最主要天氣類型。

雖然構建的天目山重污染天氣預測模型取得了較好的預報效果，但缺少要素個數、訓練與測試比率等對預測模型影響因子的量化分析，對預測模型的穩定和推廣有較大的制約。另外，研究更小尺度，如1 h、3 h等，以及非重污染天氣的預報模型和構建無縫預報體系是今后研究重點之一。