劉 宇,楚博策,高 峰,鄧 越

(1.中國電子科技集團(tuán)公司航天信息應(yīng)用技術(shù)重點(diǎn)實(shí)驗(yàn)室,河北 石家莊 050081；2.北京師范大學(xué) 地理科學(xué)學(xué)部,北京 100875)

0 引言

PM2.5是指空氣動力學(xué)當(dāng)量直徑≤2.5 μm的顆粒物，其在空氣中的濃度越高，則意味著空氣污染越嚴(yán)重。已經(jīng)有研究指出，人體某些病癥的發(fā)生概率與人在污染空氣中的暴露時(shí)間呈正相關(guān)[1-2]。隨著近年來中國北方霧霾天氣的頻繁發(fā)生，尤其是京津冀地區(qū)，社會對PM2.5的分布情況越來越關(guān)注[3]。地面PM2.5濃度的監(jiān)測及實(shí)時(shí)公布已經(jīng)成為環(huán)境領(lǐng)域的重要內(nèi)容。傳統(tǒng)的PM2.5監(jiān)測方式是地面站點(diǎn)監(jiān)測,“十二五”以來，已經(jīng)在全國建立了空氣質(zhì)量監(jiān)測網(wǎng)，但是監(jiān)測站點(diǎn)基本上集中在城市地區(qū)，且全國總共只有1 497個(gè)監(jiān)測站點(diǎn)。站點(diǎn)數(shù)量嚴(yán)重不足以及空間分布不均導(dǎo)致難以在面尺度上進(jìn)行全國PM2.5監(jiān)測[4]。

低分辨率遙感能夠獲得大范圍的大氣觀測信息，且擁有較高的時(shí)間分辨率，能夠滿足對PM2.5進(jìn)行時(shí)間和空間上的連續(xù)觀測，所以近些年已有很多研究人員使用遙感進(jìn)行不同地區(qū)的PM2.5反演[5-7]。根據(jù)PM2.5反演方法的差異，可以將這些PM2.5反演模型分為2類：基于模擬和基于衛(wèi)星觀測數(shù)據(jù)的模型?；谀M的模型通常使用全球化學(xué)傳輸模型來刻畫氣象因子和氣溶膠對PM2.5濃度的影響效應(yīng)，然后將建立的模型應(yīng)用到衛(wèi)星觀測數(shù)據(jù)，來反演PM2.5濃度[8-9]。該類模型的缺點(diǎn)是對數(shù)據(jù)要求較高，且化學(xué)傳輸模型需要較多的參數(shù)設(shè)定?；谛l(wèi)星觀測數(shù)據(jù)的模型是遙感領(lǐng)域較為通用的模型，其基本原理是基于PM2.5濃度和氣溶膠光學(xué)厚度(AOD)間的統(tǒng)計(jì)關(guān)系。該模型是遙感領(lǐng)域較為常用的方法，優(yōu)點(diǎn)是所需數(shù)據(jù)教少，且方法簡單[10]。利用地理加權(quán)回歸算法建立MODIS AOD數(shù)據(jù)與PM2.5濃度的統(tǒng)計(jì)關(guān)系，并考慮氣象因子在模型中的有效性，該模型的精度明顯優(yōu)于多元線性回歸模型[11]，Ma等[12]則同時(shí)考慮了氣象因子和土地利用情況對PM2.5反演的影響。楊麗娟等[13]認(rèn)為雖然氣象因子和土地利用情況等參數(shù)對AOD和PM2.5間的模型建立有積極影響，但氣象因子的變化較復(fù)雜，會導(dǎo)致AOD和PM2.5濃度間呈現(xiàn)明顯的日差異關(guān)系，故建立了包含固定效應(yīng)和隨機(jī)效應(yīng)的日校正模型?？紤]到PM2.5和AOD數(shù)據(jù)間具有時(shí)空變異的統(tǒng)計(jì)特性后，有研究人員嘗試使用2層統(tǒng)計(jì)模型來分別對PM2.5-AOD在時(shí)間和空間上的變異進(jìn)行建模[14-16]。半經(jīng)驗(yàn)?zāi)Ｐ褪橇硗庖环N較常用的PM2.5反演模型，該模型考慮的因子與上述模型并無差別，只是一般使用指數(shù)函數(shù)形式進(jìn)行氣象因子和PM2.5間的回歸[7,17-18]。

傳統(tǒng)的回歸模型能夠通過有限樣本刻畫簡單的映射關(guān)系，但實(shí)際應(yīng)用中，數(shù)據(jù)量的增加會給回歸模型帶來巨大挑戰(zhàn)。當(dāng)因變量與自變量間是復(fù)雜的非線性映射關(guān)系時(shí)，如PM2.5濃度與氣象因子、AOD間呈難以簡單描述的非線性關(guān)系，傳統(tǒng)的回歸模型很難對樣本進(jìn)行準(zhǔn)確建模[19]。而深度神經(jīng)網(wǎng)絡(luò)(DNN)能夠有效學(xué)習(xí)樣本中的復(fù)雜非線性映射關(guān)系，對噪聲數(shù)據(jù)有較好的魯棒性，非常適用于PM2.5反演的建模。

京津冀地區(qū)是我國霧霾發(fā)生最頻繁以及最嚴(yán)重的區(qū)域，監(jiān)測京津冀地區(qū)及其周邊PM2.5的時(shí)空分布可為環(huán)境部門提供決策支撐。當(dāng)前，使用深度學(xué)習(xí)進(jìn)行PM2.5反演建模的研究還很少；為此，本文使用MODIS AOD數(shù)據(jù)，并結(jié)合氣象同化資料，利用DNN挖掘PM2.5和氣象因子、AOD間的非線性關(guān)系，構(gòu)建PM2.5反演模型，并依此分析PM2.5的時(shí)空分布模式。

1 研究區(qū)和數(shù)據(jù)介紹

1.1 研究區(qū)概況

京津冀地區(qū)是我國空氣污染最嚴(yán)重的地區(qū)之一，根據(jù)環(huán)保部的數(shù)據(jù)顯示，2016年京津冀區(qū)域大氣優(yōu)良天數(shù)比例為56.8%，比全國平均天數(shù)比例低22%。京津冀大氣污染嚴(yán)重的原因主要有2個(gè)：其一，京津冀地區(qū)是重化工業(yè)集中的重要區(qū)域，集中了多種廢氣、廢水排放，從而直接影響大氣質(zhì)量的產(chǎn)業(yè)，其中以鋼鐵產(chǎn)業(yè)為主，包括船舶、水泥等重污染產(chǎn)業(yè)；其二，京津冀地區(qū)地形復(fù)雜，西鄰太行山脈、北靠燕山山脈、東鄰渤海，再考慮到京津冀地區(qū)的海陸風(fēng)轉(zhuǎn)換情況，京津冀地區(qū)的大氣污染物難以疏散[20]。

隨著社會對大氣污染問題的重視程度逐漸提高，國家相繼出臺了若干細(xì)則進(jìn)行京津冀地區(qū)的大氣污染整治。為及時(shí)評估大氣污染整治情況，有必要及時(shí)對京津冀地區(qū)PM2.5的時(shí)空分布進(jìn)行長期監(jiān)測。

1.2 數(shù)據(jù)介紹

本文使用的數(shù)據(jù)包括2015年中國環(huán)境監(jiān)測總站發(fā)布的PM2.5站點(diǎn)數(shù)據(jù)、MODIS AOD氣溶膠產(chǎn)品數(shù)據(jù)以及再分析氣象資料MERRA-2。

1.2.1 PM2.5站點(diǎn)數(shù)據(jù)

2013年開始，中國環(huán)境監(jiān)測總站在全國338個(gè)地級以上城市設(shè)置空氣質(zhì)量監(jiān)測站點(diǎn)1 436個(gè)，“十二五”期間，又建成了農(nóng)村區(qū)域空氣質(zhì)量監(jiān)測站點(diǎn)61個(gè)。目前，中國環(huán)境監(jiān)測總站實(shí)時(shí)公布1 497個(gè)監(jiān)測站點(diǎn)的監(jiān)測數(shù)據(jù)(包括PM2.5、PM10、NO2、SO2、O3、CO濃度)。該數(shù)據(jù)為遙感反演面尺度的PM2.5濃度提供了數(shù)據(jù)基礎(chǔ)。

由于單一省份內(nèi)空氣質(zhì)量監(jiān)測站點(diǎn)數(shù)量較少，且地理分布不均勻，會導(dǎo)致深度學(xué)習(xí)樣本過少，故本研究使用包括北京市、天津市、河北省及其周邊共9個(gè)省市的空氣質(zhì)量監(jiān)測站點(diǎn)數(shù)據(jù)，以增加樣本數(shù)量。提取2015年1月1日—12月31日的PM2.5逐小時(shí)數(shù)據(jù)，并求取10—11時(shí)的PM2.5濃度平均值，以對應(yīng)MODIS上午的過境時(shí)間10:30。

1.2.2 MODIS AOD產(chǎn)品

MOD/MYD04(MODIS Terra/Aqua Aerosol)產(chǎn)品是NASA發(fā)布的Level 2級氣溶膠產(chǎn)品，用來獲取全球海洋和陸地的大氣氣溶膠光學(xué)特性(如AOD等)。在之前的Collection 5中，NASA僅提供了10 km分辨率的氣溶膠產(chǎn)品，而在最新的Collection 6中，NASA提供了3 km分辨率的氣溶膠產(chǎn)品。本文使用2015年每天的MOD/MYD AOD產(chǎn)品，空間分辨率為3 km。

Terra星過境時(shí)間大約在10:30，Aqua星過境時(shí)間大約在13:30，故MOD04 AOD產(chǎn)品和MYD04 AOD產(chǎn)品分別對應(yīng)的是10:30和13:30的AOD數(shù)據(jù)。由于該產(chǎn)品中部分像素值缺失，故本文使用回歸方式對缺失值進(jìn)行補(bǔ)充[21]：將每天的MOD04 AOD和MYD04 AOD數(shù)據(jù)進(jìn)行線性回歸，然后利用該回歸方程進(jìn)行MOD04 AOD缺失值的補(bǔ)充。

1.2.3 氣象數(shù)據(jù)MERRA-2

本文使用NASA發(fā)布的MERRA-2再分析產(chǎn)品，該產(chǎn)品使用GEOS-5同化系統(tǒng)生成，空間分辨率為0.625°×0.5°，包括自1981年至今的氣象再分析資料。提取京津唐地區(qū)2015年每日10時(shí)的地表氣壓、2 m高相對濕度、2 m高氣溫?cái)?shù)據(jù)，行星邊界層高度數(shù)據(jù)以及地表風(fēng)速數(shù)據(jù)用于本研究的建模(GMAO,2015)。

將以上數(shù)據(jù)處理為時(shí)間和空間上統(tǒng)一的數(shù)據(jù)集。首先，將AOD數(shù)據(jù)和氣象數(shù)據(jù)重投影到地理坐標(biāo)系下，分辨率為0.03°；然后，提取對應(yīng)空氣質(zhì)量監(jiān)測站點(diǎn)位置的AOD和氣象數(shù)據(jù)，考慮到空氣質(zhì)量監(jiān)測站點(diǎn)和其他數(shù)據(jù)間的地理定位誤差，使用3×3窗口像元平均值為對應(yīng)空氣質(zhì)量監(jiān)測站點(diǎn)的AOD和氣象數(shù)據(jù)。最后，刪除無觀測值記錄，共得到32 753條記錄用來發(fā)展PM2.5反演模型，每條記錄包括PM2.5濃度、AOD、相對濕度RH、地表氣壓SP、氣溫T、行星邊界層高度PBLH及地表風(fēng)速SWS。

2 PM2.5濃度反演方法

2.1 樣本均衡化

PM2.5濃度天數(shù)呈伽馬分布如圖1所示，污染主要集中在低于100 μg/m3的區(qū)域，PM2.5濃度高的天數(shù)非常少，濃度極高的天數(shù)趨近于0。不均衡的樣本分布會導(dǎo)致DNN模型訓(xùn)練的時(shí)候難以捕獲PM2.5高值區(qū)的特征[22]，故需要對樣本進(jìn)行均衡化處理。

圖1 2015年研究區(qū)內(nèi)所有空氣質(zhì)量監(jiān)測站點(diǎn)每日10時(shí)PM2.5濃度天數(shù)分布Fig.1 Daily distribution of PM2.5 concentrations at 10 a.m.observed by all air quality monitoring systems in the study region in 2015

本文使用傳統(tǒng)的重復(fù)過采樣方法進(jìn)行樣本均衡化。假設(shè)樣本集中PM值為i的天數(shù)為Ni，首先確定Ni的最大值Nmax，然后將Ni小于Nmax的樣本復(fù)制至Nmax個(gè)，使PM濃度值呈均勻分布。

2.2 深度神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)及其評估方法

將均衡化的樣本按照0.9∶0.1的比例劃分為訓(xùn)練樣本與檢驗(yàn)樣本；根據(jù)訓(xùn)練樣本與檢驗(yàn)樣本不斷調(diào)整DNN的結(jié)構(gòu)參數(shù)，包括隱藏層數(shù)量和每層的節(jié)點(diǎn)數(shù)量，得到的最優(yōu)結(jié)構(gòu)如圖2所示。該DNN共有8個(gè)隱藏層，每層的隱藏節(jié)點(diǎn)個(gè)數(shù)分別為39,31,27,21,15,11,7,3。

圖2 擬合效果最優(yōu)的深度神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.2 Structure of DNN obtaining optimal fitting result

作為對比，使用式(1)的半經(jīng)驗(yàn)回歸模型對樣本進(jìn)行回歸驗(yàn)[18]：

ln(PM2.5)=β0+βAODln(AOD)+βPBLHln(PBLH)+

βRHln(RH)+βTln(T)。

(1)

對該模型的評估指標(biāo)包括相關(guān)系數(shù)r(PM2.5真實(shí)值與模型值間的相關(guān)系數(shù))、平均絕對誤差MAE以及平均相對誤差MRE。

2.3 基于點(diǎn)面融合模型的PM2.5時(shí)空數(shù)據(jù)融合方法

由于AOD數(shù)據(jù)缺失，導(dǎo)致DNN反演得到的PM2.5分布也是缺失的，本文采用點(diǎn)面融合模型對PM2.5的時(shí)空分布進(jìn)行插值[21]?？諝赓|(zhì)量監(jiān)測站點(diǎn)的PM2.5濃度數(shù)據(jù)不包含空間特征，但其反映的時(shí)序規(guī)律是準(zhǔn)確的，而DNN反演的PM2.5濃度數(shù)據(jù)則包含部分準(zhǔn)確的空間特征；點(diǎn)面融合模型結(jié)合了站點(diǎn)插值數(shù)據(jù)與DNN反演數(shù)據(jù)的優(yōu)點(diǎn)，該模型假設(shè)站點(diǎn)插值PM2.5數(shù)據(jù)在不同日期間的差值與DNN反演的PM2.5數(shù)據(jù)在相同日期間的差值相同，具體方法為：

Vi,j,t1-Vi,j,t2=Mi,j,t1-Mi,j,t2，

(2)

式中，i,j代表像素位置(行、列值)；t代表時(shí)間(天)；V為DNN反演的PM2.5濃度數(shù)據(jù)；M為PM2.5站點(diǎn)插值數(shù)據(jù)。假設(shè)Vi,j,t1為PM2.5濃度缺失值，Vi,j,t2為時(shí)間上距離Vi,j,t1最近的有效值，則其差值與PM2.5濃度在站點(diǎn)插值數(shù)據(jù)中的差是相同的。

3 結(jié)果與討論

3.1 樣本均衡化的效果

原始樣本數(shù)量為32 753個(gè)，經(jīng)過樣本均衡化后的樣本數(shù)量為148 330個(gè)。對原始樣本同樣按照0.9∶0.1的比例進(jìn)行DNN的訓(xùn)練和檢驗(yàn)，得到擬合效果最優(yōu)的網(wǎng)絡(luò)結(jié)構(gòu)與圖2相同。作為均衡化樣本的對比，表1列出了原始樣本與均衡化樣本的擬合結(jié)果。

表1 DNN模型對原始樣本與均衡化樣本擬合效果對比Tab.1 Comparison between fitting result for original samples and balanced samples separately on DNN model

原始樣本的擬合相關(guān)系數(shù)0.43遠(yuǎn)低于均衡化樣本的擬合精度0.94；同時(shí)由于原始樣本分布過于集中，導(dǎo)致訓(xùn)練樣本集與檢驗(yàn)樣本集的擬合效果相差較大，訓(xùn)練樣本集的擬合相關(guān)系數(shù)為0.68，而檢驗(yàn)樣本集的擬合相關(guān)系數(shù)僅為0.43。原始樣本與均衡化樣本擬合結(jié)果的MAE差異不大，但由于原始樣本分布過于集中，其訓(xùn)練樣本集與檢驗(yàn)樣本集的MAE差異更大；原始樣本中訓(xùn)練樣本集與檢驗(yàn)樣本集擬合結(jié)果的MRE差異巨大，達(dá)到20.6%。

作為對比，均衡化樣本的擬合結(jié)果遠(yuǎn)遠(yuǎn)優(yōu)于原始樣本，同時(shí)訓(xùn)練樣本和檢驗(yàn)樣本在模型擬合精度上差異也明顯優(yōu)于原始樣本，MAE差異僅為0.5 μg/m3，MRE差異僅為0.3%，故樣本均衡化明顯提高了DNN模型的魯棒性和PM2.5反演時(shí)的精度。

3.2 PM2.5回歸精度

使用圖2的DNN模型分別對訓(xùn)練樣本集和檢驗(yàn)樣本集進(jìn)行擬合，擬合結(jié)果如圖3所示。對訓(xùn)練樣本集而言，真實(shí)PM2.5濃度與模型預(yù)測的PM2.5濃度間的相關(guān)系數(shù)為0.944，平均絕對偏差為24.5 μg/m3，平均相對偏差為25.9%；對于檢驗(yàn)樣本集而言，真實(shí)PM2.5濃度與模型預(yù)測的PM2.5濃度間的相關(guān)系數(shù)為0.94，平均絕對偏差為25.0 μg/m3，平均相對偏差為26.2%，DNN對訓(xùn)練樣本集和檢驗(yàn)樣本集的回歸散點(diǎn)圖如圖3所示。

(a) 訓(xùn)練樣本集的擬合結(jié)果

半經(jīng)驗(yàn)回歸模型對訓(xùn)練樣本集和驗(yàn)證樣本集的集合結(jié)果如表2所示，回歸相關(guān)系數(shù)僅為0.18，遠(yuǎn)低于DNN模型的0.94；平均絕對誤差和平均相對誤差也遠(yuǎn)高于DNN模型的結(jié)果。

表2 半經(jīng)驗(yàn)回歸模型對原始樣本與均衡化樣本擬合效果對比Tab.2 Comparison between fitting result for original samples and balanced samples separately onsemi-empirical model

3.3 PM2.5時(shí)空分布

利用點(diǎn)面融合算法，可以得到每天上午10:30空間分辨率為0.03°的PM2.5時(shí)空分布；考慮到在氣象條件不發(fā)生較大變化的情況下，PM2.5濃度不會發(fā)生劇烈變化，故上午10:30的PM2.5濃度可以一定程度上代表全天的PM2.5濃度。本文假定每日10:30的PM2.5濃度為日均PM2.5濃度，全月每日10:30的平均PM2.5濃度為月均PM2.5濃度，全年每日10:30的平均PM2.5濃度為年均PM2.5濃度。

DNN反演得到的月均PM2.5時(shí)空分布如圖4所示。

圖4 DNN模型反演的京津冀地區(qū)2015年月均PM2.5濃度時(shí)空分布圖Fig.4 Spatiotemporal distribution of average PM2.5 concentration monthly in Beijing-Tianjin-Hebei area retrieved by DNN model

從時(shí)間上看，PM2.5濃度較高的月份主要集中在屬于冬季的12月、1月和2月，這與實(shí)際情況相符，每年冬季由于京津冀地區(qū)供暖的需求，煤炭燃燒會向大氣輸送大量污染物，同時(shí)，受蒙古西伯利亞高壓控制，大氣污染物不易擴(kuò)散；PM2.5濃度最低的月份為9月。

PM2.5濃度年均地理分布如圖5所示。PM2.5濃度較高的地區(qū)集中在北京市中部以南地區(qū)，太行山東部，與實(shí)際情況相符，受到地形因素影響，京津冀地區(qū)西側(cè)毗鄰太行山，北京北部山地環(huán)繞，導(dǎo)致大氣污染物難以擴(kuò)散；京津冀北部空氣質(zhì)量為優(yōu)良狀況。

(a) 空氣質(zhì)量監(jiān)測站點(diǎn)插值數(shù)據(jù) (b) DNN模型反演數(shù)據(jù)圖5 PM2.5年平均空間分布Fig.5 Annual average geographical distribution of PM2.5 concentration

作為京津冀地區(qū)空氣污染嚴(yán)重且最受關(guān)注的城市，北京市與石家莊市PM2.5濃度的逐日波動如圖6所示。

(a) 北京市

對于北京市，DNN反演的11、12月份PM2.5濃度與空氣質(zhì)量監(jiān)測站點(diǎn)數(shù)據(jù)一致性很好，但其他日期的數(shù)據(jù)間存在系統(tǒng)偏差；對于石家莊市，DNN反演的全年P(guān)M2.5濃度與空氣質(zhì)量監(jiān)測站點(diǎn)間數(shù)據(jù)一致性都很好。北京市全年的空氣污染波動頻率較均勻，PM2.5濃度峰值主要分布在冬季，最高值可達(dá)589 μg/m3；12月份的空氣污染頻率略低于其他月份，但大氣污染的持續(xù)時(shí)間長于其他月份，年平均PM2.5濃度為107.3 μg/m3。石家莊市呈現(xiàn)出與北京市相同的空氣污染波動模式，PM2.5濃度峰值同樣主要分布在冬季，但峰值395 μg/m3低于北京市，但年平均PM2.5濃度110.5 μg/m3略高于北京市；同時(shí)，石家莊市冬季高污染的頻率也高于北京市。

4 結(jié)束語

本文使用MODIS AOD數(shù)據(jù)及MERRA2氣象再分析資料構(gòu)建了PM2.5濃度反演的DNN模型，驗(yàn)證了該模型在PM2.5濃度反演方面的效果，通過樣本均衡化的方法提高了DNN在反演PM2.5濃度時(shí)的精度，并用此模型反演了2015年京津冀地區(qū)的PM2.5濃度時(shí)空分布。最后，通過月均PM2.5濃度空間分布、年均PM2.5濃度空間分布及北京市與石家莊市PM2.5日波動情況進(jìn)行了初步分析。得到相關(guān)結(jié)論：① 樣本均衡化增加了PM2.5高值的樣本數(shù)量，使得樣本分布更加均一，有效地將DNN的擬合相關(guān)系數(shù)由0.43提高到了0.94，且增加了DNN模型的魯棒性。DNN模型對PM2.5濃度站點(diǎn)數(shù)據(jù)的擬合相關(guān)系數(shù)分別為0.944和0.94，平均絕對誤差和平均相對誤差分別為25.0 μg/m3，26.2%，完全能滿足對京津冀地區(qū)大氣污染監(jiān)測的需要。② 京津冀地區(qū)的大氣污染主要分布在中南部，且中南部的污染程度遠(yuǎn)高于北部地區(qū)；2015年，京津冀地區(qū)大氣污染最嚴(yán)重的時(shí)間為冬季，夏季的PM2.5濃度最低，北京市與石家莊市PM2.5濃度峰值均分布在冬季；北京市與石家莊市12月份的空氣污染頻率略低于其他月份，但大氣污染的持續(xù)時(shí)間長于其他月份。