基于混合效應模型的京津冀地區PM2.5數值模擬*

鮑艷磊 田 冰# 胡引翠 張文靜 王晨旭 王天宇

(1.河北師范大學資源與環境科學學院，河北 石家莊 050024；2.河北省環境演變與生態建設實驗室，河北 石家莊 050024)

自20世紀90年代以來，霧霾逐漸成為中國嚴重的環境問題之一[1]。PM2.5是霧霾的主要污染物，許多流行病學研究表明，PM2.5對人體健康產生不良影響，造成呼吸、心血管、免疫和生殖系統出現功能障礙，嚴重甚至會導致死亡[2-4]。近年來，其引起了國家和社會的廣泛關注。目前，PM2.5自動監測依賴于地面監測格網，由于目前中國監測格網數量還有所欠缺，對于研究區域的PM2.5濃度獲取存在困難，而利用遙感衛星獲得的氣溶膠光學厚度(AOD)數據反映近地面污染情況具有連續性好、傳輸速度快等難以比擬的優勢。以AOD為介質的消光系數在垂直方向上的積分，是表征大氣渾濁程度的關鍵物理量，也是確定氣溶膠氣候效應的重要因素。通常，高的AOD預示著氣溶膠縱向積累的增長，直接導致大氣能見度的降低。國內外許多學者開展了利用AOD估算PM2.5的研究。研究表明，AOD和PM2.5之間存在明顯的正相關關系。近年來，不同學者基于多元線性模型[5]、土地利用回歸模型[6-9]、協同時空地理加權回歸模型[10]等估算PM2.5，基本得到了很好的結果。由此可知，利用AOD來估算地面PM2.5是可行的，對于了解PM2.5的時空分布及彌補沒有監測站點地區的數據缺失是具有參考價值的。

2018年，北京市、天津市、河北全省PM2.5年平均值分別約為51、52、56 μg/m3，京津冀地區平均值超過了《環境空氣質量標準》(GB 3095—2012)一級標準(35 μg/m3)。針對該地區大氣質量嚴重污染問題，《大氣污染防治行動計劃》提出建立京津冀等三大區域大氣污染聯防聯控機制，開啟京津冀區域協作治理的新進程。

郝靜等[11]使用空間分辨率為10 km的AOD數據，利用混合效應模型對京津冀內陸平原PM2.5濃度進行時空變化定量模擬，得到R2=0.78。景悅等[12]2890利用基于每日的混合效應模型探究AOD與PM2.5的相關關系，R2=0.92。由此說明，高精度的遙感數據和混合效應模型為預測PM2.5濃度提供了更好的數據源和數據模型。但PM2.5與AOD的關系會受到氣象因素和其他因素的影響而存在不穩定性。國內外眾多學者利用AOD來估算PM2.5濃度時加入氣象要素，結果發現，這些因素會在一定程度提高兩者相關性[13-15]。薛文博等[16]研究表明，PM2.5與訂正過的AOD的R2為0.77；同時引入相對濕度、風速，近地面PM2.5與反演值的R2升至0.79。目前，在混合效應模型中引入氣象因子作為自變量對京津冀地區PM2.5與AOD相關關系的研究較少，本研究嘗試以AOD數據和氣象要素為預測因子，得到最優的混合效應模型，以期達到更好預測京津冀地區PM2.5濃度的效果。

1 研究區概況

京津冀是中國的“首都經濟圈”，包括北京市、天津市及河北省11個地級市。京津冀地區包括多種類型的地貌，包括西部的壩上高原、燕山、太行山脈和中部、東南部的平原地區。京津冀地區屬于溫帶季風性氣候，具有明顯的季節變化特征。燃煤、機動車、工業廢氣排放多種污染物共生局面加之西、北面環山的地形不利于霧霾的擴散情況導致近年來京津冀地區大氣污染超標頻度較高。為更好監測PM2.5對大氣環境影響，目前京津冀地區有PM2.5監測站點80個、氣象站點26個，具體位置如圖1所示。

2 數據與方法

2.1 數據來源

2.1.1 中分辨率成像光譜儀(MODIS)AOD數據

MODIS數據通過將MODIS搭載在Terra衛星上獲得，光譜范圍0.4～14.4 μm。本研究采用美國國家航空航天局(NASA)網站(https：//ladsweb.nascom.nasa.gov)發布的以日為時間分辨率、3 km為空間分辨率的MOD04 L2 C6產品數據，時間為2017年1月1日至12月31日。

首先在ENVI5.1軟件中對下載的原始數據進行幾何校正、拼接和裁剪等預處理，之后在ArcGIS 10.2軟件中獲取PM2.5監測站點周圍3×3個柵格像元的AOD平均值與對應PM2.5監測站點相匹配。

2.1.2 地面PM2.5數據

京津冀地區的80個PM2.5監測站點的數據來自中國環境監測中心官方網站(http：//106.37.208.233：20035/)。Terra衛星在中國的過境時間是10：30，因此本研究采用衛星過境時段(即10：00—11：00)地面監測站點的PM2.5平均值。在以3 km為單元劃分的京津冀地區監測格網中可能會存在一個或多個PM2.5監測站點，此時PM2.5取監測格網中所有監測站點的平均值。

2.1.3 氣象和邊界數據

京津冀地區26個氣象站點的氣象數據可從中國氣象科學數據共享服務網(http：//data.cma.cn/)獲得，包括相對濕度、風速、氣溫和氣壓；因氣象站點與PM2.5監測站點經緯度不同，因此需要將氣象站點的數據進行插值得到與PM2.5監測站點相匹配的氣象數據。

省、地、市行政邊界數據來源于國家基礎地理信息系統網站(http：//ngcc.sbsm.gov.cn/)。

2.2 研究方法

2.2.1 混合效應模型

混合效應模型是指既包含固定效應又包括隨機效應的模型。固定效應類似于標準回歸系數，直接估計得到。隨機效應從它們的方差和協方差估計值中總結而來。隨機效應用于解釋AOD數據、氣象因子與PM2.5之間的日變化關系，以隨機截距或隨機系數的形式呈現。

景悅等[12]2892-2894指出，基于日變化的混合效應模型能在時間上解釋PM2.5和AOD的關系變化，基于站點的混合效應模型也能在一定程度上解釋空間上PM2.5和AOD的關系變化，而且兩種模型均優于簡單的線性回歸。式(1)表示在自變量只有AOD的情況下建立的3種模型：(1)站點變化影響下的混合效應模型(基于空間的模型)；(2)時間變化影響下的混合效應模型(基于時間的模型)；(3)PM2.5和AOD關系在時間和站點共同影響下的混合效應模型(基于時空的模型)。3種模型使用相同的數據，基于不同的變量字段來實現。式(2)表示增添了氣象因子作為影響因子的3種混合效應模型，反映PM2.5和AOD關系在不同情況下的變化。

圖1 京津冀PM2.5監測站點和氣象站點分布Fig.1 Distribution of Beijing-Tianjin-Hebei PM2.5 monitoring site and meteorological site

Mij=(α+uj)+(β1+νj)×Aij+Si+εij

(1)

Mij=(α+uj)+(β1+νj)×Aij+(β2+κj)×Tij+(β3+lj)×Wij+(β4+mj)×Rij+β5×Pij+Si+εij

(2)

式中：Mij為第i個監測格網第j天的PM2.5，μg/m3；α為固定截距；uj為第j天隨機截距；β1、β2、β3、β4、β5為與AOD、氣溫、風速、相對濕度、氣壓有關的固定斜率；vj、κj、lj、mj為第j天與AOD、氣溫、風速、相對濕度、氣壓有關的隨機效應斜率；Aij為第i個監測格網第j天的AOD；Si為第i個監測格網的隨機效應；εij為第i個監測格網第j天的隨機誤差項；Tij為第i個監測格網第j天的氣溫，℃；Wij為第i個監測格網第j天的風速，m/s；Rij為第i個監測格網第j天的相對濕度，%；Pij為第i個監測格網第j天的氣壓，Pa。

2.2.2 模型驗證方法

采用十折交叉驗證法驗證混合效應模型的準確度。具體實現過程：將數據集隨機分成10等份，輪流將9份參與建模，剩余一份作為測試集，重復10次直至所有數據都被預測出來，然后將10次驗證結果的均值作為對預測模型精度的估計。

模型擬合效果的判斷以PM2.5監測值與模擬值的R2、十折交叉驗證后相關系數(CVR2)、均方根誤差(RMSE，μg/m3)及平均絕對誤差(MAE，μg/m3)表示。

3 結果與討論

3.1 模型擬合

采用基于不同類型的混合效應模型，并嘗試在固定效應和隨機效應中增添、改變氣象因子的類型和數量，來進行模型的擬合和驗證，共嘗試進行了78種類型的含有不同自變量的混合效應模型。

在固定效應和隨機效應中加入不同氣象因子后，R2、CVR2、RMSE、MAE的取值見表1。R2和CVR2取值越大，兩者的相關性越強；RMSE和MAE取值越小，模型的誤差越小。綜合分析，基于空間的模型的R2、CVR2比另外兩種模型小，基于時間的模型的R2、CVR2基本都比基于時空的模型小；基于空間的模型的RMSE、MAE明顯大于另外兩種模型，后兩種模型數值有交叉，但總體基于時空的模型的RMSE、MAE小于基于時間的模型。

在78種模型中，分別選擇基于空間、時間和時空的最優模型進行分析。基于空間的模型，在固定效應和隨機效應中同時含有AOD、風速和相對濕度的模型得到的結果最優；基于時間的模型，在固定效應和隨機效應中同時含有AOD、風速、相對濕度、氣溫和氣壓的模型的結果最優；基于時空模型，在固定效應和隨機效應中同時含有AOD、風速、相對濕度和氣溫，且僅在固定效應中含有氣壓的模型的結果最優。將這3種最優模型經十折交叉模型校正，結果見圖2。模型R2=0.90，校正后基于時空最優模型擬合程度最好，CVR2=0.81，RMSE、MAE分別為13.44、10.12 μg/m3。表1中，基于時空的模型的CVR2比R2小，說明模型存在部分過度擬合，但總體來說模型的顯著性明顯，說明氣象因子可提高模型的擬合精度和PM2.5的預測精度。

表1 R2、CVR2、RMSE、MAE的取值

圖2 模型模擬值與監測值間校正后得到的相關關系Fig.2 Correlation diagram between predicted and actual values of the model

基于時空的模型中引入氣溫、相對濕度、風速和氣壓4種氣象因子，此模型得到的固定效應(包括固定截距和固定斜率)的解見表2。其中，混合效應模型的固定效應類似于標準回歸系數，表示全年內自變量對于PM2.5的固定影響，AOD固定斜率估計值為21.861 6，說明PM2.5與AOD呈現正相關；相對濕度固定斜率估計值0.371 0，說明相對濕度越大越有利于大氣顆粒物質的形成，PM2.5濃度越大；氣溫固定斜率斜率估計值為0.064 2，說明PM2.5與氣溫呈現正相關，氣溫通過控制邊界層高度影響大氣顆粒物質的擴散；氣壓固定斜率估計值為0.000 1，說明氣壓增大，PM2.5濃度也提高，但影響程度較小。風速固定斜率估計值為-0.453 2，表明PM2.5濃度與風速呈反比，風有利于大氣顆粒物的擴散，風速越大越有利于空氣中的PM2.5的輸送與擴散。風速與相對濕度的系數較高，是相對主要的氣象影響因子。

由表3可知，各城市的CVR2為0.71～0.88，RMSE為9.79～14.92 μg/m3，MAE為8.19～11.80 μg/m3。根據CVR2，北京市、邢臺市、秦皇島市、張家口市、廊坊市和保定市的擬合效果較好，石家莊市、承德市的擬合效果較差，但總體能在一定程度上表達PM2.5與AOD的關系，表現為較高的擬合效果、較低的誤差。

3.2 京津冀地區PM2.5濃度模型適應性與時空分布分析

2017年，京津冀地區PM2.5監測值為1～232 μg/m3，平均為45.23 μg/m3；PM2.5模擬值為0～229 μg/m3，平均為44.60 μg/m3。分別將各監測站點的PM2.5監測、模擬年均值插值到整個研究區，結果見圖3。PM2.5濃度空間差異性和數值吻合度均較高，說明模型的模擬效果較好，因此可利用該模型對缺乏監測數據地區以及大范圍研究區進行PM2.5的模擬與預測。

為方便分析，本研究將京津冀分為4個地區，分別為北京市、天津市、河北平原區(廊坊市、保定市、石家莊市、衡水市、邢臺市和邯鄲市)和河北非平原區(唐山市、秦皇島市、滄州市、張家口市和承德市)。2017年，北京市、天津市、河北非平原區、河北平原區的PM2.5模擬年均值分別為42.51、51.38、38.32、46.72 μg/m3，整體上北部低于南部、東部沿海低于西部地區。

本研究還將一年分為冷暖兩季(冷季為10月至次年3月，暖季為4—9月)，對PM2.5濃度及模型適應性進行分析，結果見圖4。京津冀全區冷、暖季PM2.5模擬平均值分別為44.99、44.49 μg/m3。全區暖、冷季PM2.5模擬值與監測值分別相差0.31、0.36 μg/m3，其中北京市分別相差0、1.76 μg/m3，天津市分別相差0.17、2.00 μg/m3，河北平原區分別相差0.35、1.59 μg/m3，河北非平原區分別相差0.83、1.10 μg/m3。全區平均值相對誤差值不超過0.008，分區平均值相對誤差值為0～0.037。總體來說，暖季相較于冷季估計均值更精確，模型對PM2.5濃度估計精度較高。

表2 固定效應的解

表3 基于城市統計的混合效應模型擬合結果及校正結果

圖3 2017年PM2.5監測與模擬年均值分布對比Fig.3 Comparison of the PM2.5 distribution of monitoring values and simulation values in 2017

圖4 基于不同地區的PM2.5監測值與模擬值對比Fig.4 Comparison of PM2.5 monitoring values and simulated values in different regions

研究區4個地區PM2.5模擬年均值均大于35 μg/m3，但是不同地區也存在差別。河北非平原區中的張家口市、承德市PM2.5年均值較低，分析原因可能是因為當地工廠相對較少，空氣流通狀況好；秦皇島市為沿海城市，常年有海風的作用，因此PM2.5年均值相對較低，空氣質量較好。河北平原區工廠多、人口密集，產生PM2.5的源頭多；冷季受西北風影響，隨風飄蕩的污染物在遇到太行山時無法翻越，造成累積，然后便有可能形成氣團污染到其他周邊城市，造成二次及三次污染，導致內部平原地區PM2.5年均值明顯偏高。北京市人口密集，但工廠外遷，冬季盛行西北風，濕度較小，較有利于PM2.5的擴散，所以PM2.5年均值相對不是很高。整體來說，京津冀地區呈現南高北低、平原地區東南高西北低、非平原地區高緯度區域相對低的空間格局，京津冀地區空氣質量亟待解決。

4 結 語

(1) MODIS 3 km AOD的數據全球開放、獲取方便，因此本研究使用此數據來預測地面的PM2.5濃度。目前，3 km的AOD產品可通過MODIS的暗像元算法(DT)和深藍算法(DB)得到，但通過DB得到的產品沒有發布，應用DT在冷季云、霧及重污染天氣反演得到的AOD產品質量較差，導致部分冬季天氣惡劣地區參與建模的數據較少。在以后的研究中可嘗試將3、10 km數據進行融合或采用SARA算法反演出500 m分辨率的AOD產品填補AOD的缺失，以提高研究的精度。另外，以后的研究中還可加入縣級監測站點，增加站點數量與均勻程度進一步提高模型的精度。

(2) 局地氣象是影響PM2.5濃度的重要因子，每天每個地方的氣象因子存在差異，因此可通過最優插值得到監測站點的各種氣象因子的數值。在構建模型時加入不同的氣象因子，所得的CVR2存在不同。其中，風速和相對濕度的影響系數較大。

(3) 基于時空的模型擬合結果比基于空間、時間的模型優。

(4) 整體來說，京津冀地區呈現南高北低、平原地區東南高西北低、非平原地區高緯度區域相對低的空間格局，京津冀地區空氣質量亟待解決。