基于深度學習的日間逐小時地表PM2.5濃度反演

徐成康，陳斯婕，董長哲，徐文韜，劉 東

（1.浙江大學 光電科學與工程學院，浙江 杭州 310027；2.上海衛(wèi)星工程研究所，上海 200240）

0 引言

PM2.5即空氣動力學直徑小于2.5 μm 的細顆粒物，又被稱為可入肺顆粒物，是空氣污染霧霾危害中的主要成分之一，對氣候、環(huán)境和人類健康有極大的危害［1］。目前，我國已形成系統(tǒng)化的地面空氣質量監(jiān)測體系，能夠對PM2.5、PM10、一氧化碳、二氧化硫、二氧化氮和臭氧等主要污染物進行實時精確監(jiān)控。但地面監(jiān)測站點缺乏區(qū)域覆蓋能力，且分布不均，無法對各類污染物的分布及相互作用過程進行有效的追蹤分析和預測，而星載遙感觀測具有全天時、高覆蓋等獨特優(yōu)勢，在大氣狀態(tài)監(jiān)測領域中被廣泛應用。如美國航空航天局（National Aeronautics and Space Administration，NASA）的TERRA 和AQUA 兩顆衛(wèi)星上的中分辨率成像光譜儀（Moderate Resolution Imaging Spectroradiometer，MODIS），在全球氣溶膠觀測上取得了巨大成功［2］。

目前已經有相關研究將星載傳感器反演的氣溶膠光學厚度（Aerosol Optical Depth，AOD）通過經驗或半經驗統(tǒng)計、地理加權回歸、化學傳輸模式算法等方法反演地表PM2.5濃度［3-5］，但僅通過被動遙感數據反演PM2.5精度較低。星載被動遙感手段能夠進行超大范圍顆粒物柱總量探測，但一定程度上受限于其垂直探測能力及光照條件，導致觀測結果較粗糙，而地面空氣質量監(jiān)測體系具有全天時高效探測近地表顆粒物濃度的能力［6］。因此，本文通過神經網絡結合星載遙感數據、地表測量數據和其他輔助數據來反演估計PM2.5，能夠提高觀測范圍和預測能力，為相關科研目標及政策實施提供便利。隨著衛(wèi)星平臺和遙感觀測儀器的高速發(fā)展，使在衛(wèi)星上搭載的先進遙感儀器進行宏觀大氣狀態(tài)觀測和源匯分析是必然的趨勢。

1 研究區(qū)域與數據來源

1.1 研究區(qū)域

研究區(qū)域為長三角地區(qū)，以杭州市國家基準氣象站為中心，地理位置如圖1 所示。長三角地區(qū)遠離中國西北部的天然沙塵源，地勢平坦植被多，基本不存在沙塵長距離傳輸［7］。從地理上看，長三角（30°～33° N，119°～122° E）是長江與錢塘江形成的天然沖積平原，工農業(yè)十分發(fā)達，且人口稠密，上海、南京、杭州等特大城市均位于該區(qū)域。圖1 中標注了杭州市國家基準氣象站的位置，本研究使用的霧霾層高度（Haze Layer Height，HLH）由該站點的地基微脈沖激光雷達（Micro-Pulse Lidar，MPL）數據計算得到。圖1 同樣標注了該區(qū)域內空氣質量指數（Air Quality Index，AQI）站點的位置分布，是PM2.5實測數據的數據來源。在研究時間上，覆蓋了從2017年1月—2020年5月共3.5a的時間。

圖1 研究區(qū)域及MPL/AQI 站點分布Fig.1 Studied area and the locations of the MPL/AQI stations

1.2 地基數據

本文使用的PM2.5濃度數據來自于AQI 站點，在研究區(qū)域內共計有182 個站點，所采集的數據包括站點每小時及24 h 滑動平均的PM2.5濃度數據，HLH 數據從杭州市國家基準氣象站點數據中計算得到。在將AOD 通過經驗公式轉化為PM2.5濃度的過程中，混合層高度（Mixing Layer Height，MLH）作為轉換中的比例因子是一個關鍵參數。AOD 是垂直方向的氣溶膠光學總量，MLH 是近地面的氣溶膠混合高度，在沒有上層傳輸的情況下（且絕大部分為細顆粒物），可以認為AOD 和PM2.5成正比，HLH 的高度和PM2.5濃度成反比［8］。在理想狀態(tài)下，邊界層高度（Boundary Layer Height，BLH）可視為混合層高度，但實際中氣溶膠并不是總局限在邊界層以內，且在邊界層內分布也不均勻，因此使用BLH 替代MLH，將柱AOD 轉化為地表PM2.5會有顯著的誤差，需要合適的校正方法來提高PM2.5濃度估算的準確性［9-10］。已有相關研究使用梯度法從微脈沖激光雷達廓線中計算BLH 并轉換為HLH，該方法充分考慮了邊界層以上的氣溶膠，能夠有效校正偏差，提高反演準確性［11-13］。

利用地基站點的監(jiān)測數據來驗證AOD-PM2.5的估算準確率，并使用HLH 標準化來代替?zhèn)鹘y(tǒng)的AOD/BLH 標準化，進一步改善估算的結果。如圖2（a）所示，無法直接從AOD 與PM2.5濃度變化中提取兩者的關系，因此通過計算皮爾遜相關系數平方（R2）來計算特征相關性，進一步分析HLH 的相關性及相較于BLH 的進步。比較圖2（b）、圖2（c），分別對應AOD/BLH 與PM2.5的濃度相關性及AOD/HLH 與PM2.5的濃度相關性，同時考慮相對濕度（Relative Humidity，RH）對PM2.5濃度探測的影響，通過計算吸濕增長因子f（RH）進行PM2.5校正，校正后兩者的R2分別為0.17 和0.41，明顯提升。

圖2 AOD 和PM2.5時序圖及相關性比較Fig.2 Time-series diagram and correlation comparison of AOD/PM2.5

1.3 AOD 與氣象數據

遙感AOD 數據來自NASA 的MODIS 傳感器公開數據。MODIS 是搭載在TERRA 和AQUA兩顆衛(wèi)星上的MODIS，覆蓋了0.4～14.4 μm 光譜范圍，分辨率最高達250 m［14］。本研究使用的數據產品為MODIS 多角度大氣校正算法（Multi-Angle Implementation of Atmospheric Correction，MAIAC）產品，編號MCD19A2。MAIAC 產品為MODIS 觀測結果通過算法反演得到的1 km 空間分辨率的AOD 數據，使用了改進的光譜表面反射算法來確保AOD 反演的準確性［15］。

由于星載AOD 數據受到云和地表特征的影響，存在較多無效值，需要模型AOD 數據輔助，并通過神經網絡對MAIAC AOD 數據進行填補。采用的AOD 模型數據來自于為戈達德地球觀測系統(tǒng)（Goddard Earth Observing System Model，GEOS）模型，數據產品為GEOS5 FP 2 d 時間平均初級氣溶膠診斷（Time-Averaged Primary Aerosol Diagnostics）空間分辨率為0.25°×0.312 5°，時間分辨率為3 h［16］。其他氣象數據來自于歐洲中期天氣預報中心（The European Centre for Medium-Range Weather Forecasts，ECMWF）的再分析數據（ECMWF Reanalysis v5，ERA5），該數據庫提供了從1979 年起，大氣、陸地和海洋氣候變量的每小時估計值［17］。本文主要使用了ERA 5 數據庫中的單層和氣壓分層數據中的地表溫度、濕度、壓強、風速和邊界層高度數據，其空間分辨率為0.25°×0.25°，時間分辨率為1 h。

2 研究方法與模型

2.1 機器學習模型

支持向量機（Support Vector Machine，SVM）常用于處理機器學習中的分類問題，對于分類問題，SVR 是在特征空間中尋找一個曲面，使得不同特征之間的間隙最大，而離該曲面最近的特征向量，則被稱為支持向量。支持向量回歸（Support Vector Regression，SVR）是SVM 方法的擴展，簡單的線性回歸是要找出一條殘差最小的線，而SVR 則是找出一個超平面，使得所有的樣本點離超平面的總偏差最小，主要用于處理回歸問題［18］。

隨機森林（Random Forest，RF）是典型的多模型聯(lián)合的機器學習模型，其通過組合一系列決策樹的結果來進行聯(lián)合預測。具體算法流程為：先對原始訓練數據集進行隨機采樣，用于訓練第一個決策樹，之后進行重新采樣訓練下一個模型，重復訓練N個模型，采樣使用隨機子空間法，確保決策樹之間的低相關性；訓練完成后，輸入測試數據能夠得到N個預測輸出，通過多數投票或者求取平均值的方法得到最終的預測結果［19］。

多層感知機（Multilayer Perceptron，MLP）是一種神經網絡模型，由輸入、隱藏層和輸出層構成，層與層之間全連接，即每一個節(jié)點都與其上下層的所有節(jié)點有連接，每一層的輸出傳入到下一層直到生成最終的結果。MLP 使用如Tanh 或Sigmoid 等激活函數，通過給神經元引入非線性因素來克服線性模型的限制，增強網絡的適用性［20］，在分類與回歸問題中有廣泛應用。

2.2 ADRN 與集成學習

由于MAIAC AOD 觀測結果存在大量缺失現象，很難使用基于卷積神經網絡的深度學習方法，因為卷積神經網絡需要完整的圖像或具有有限隨機缺失值的圖像進行訓練，而MAIAC AOD 則通常是因為云層以及地理條件等原因產生缺失，這種缺失過程并不隨機。本文使用的深度學習網絡為基于自動編碼器的深度殘差網絡（Auto-encoderbased Deep Residual Network，ADRN）及集成算法，算法結構及流程如圖3 所示。ADRN 具有對稱結構的編碼和解碼層，能夠對多個協(xié)變量特征進行監(jiān)督學習，通過中間層潛在地學習有效的數據特征，提高學習效率［21］。普通的前向反饋神經網絡（Feedforward Neural Network，FFNN）在用于回歸預測時，會隨著隱藏層的增加而發(fā)生梯度飽和以及精度下降的現象，因此ADRN 在編碼層和解碼層引入殘差連接以提高學習效率［22］，實現在廣區(qū)域內穩(wěn)健估算AOD 數據。

圖3 ADRN 模型結構及集成學習流程Fig.3 ADRN model structure and the integrated learning process

模型訓練時使用了集成裝袋（Bagging）算法，該算法被廣泛應用于回歸和分類應用中，上文提到的RF 算法就是一種典型的集成算法。由于使用的數據集涵蓋地時間跨度較大（3.5 a），且數據分布離散，訓練單個模型受數據集的影響較大，因此預測效果較差。而Bagging 算法通過將原始訓練集進行隨機取樣，構建并訓練多個獨立的弱學習器，之后將弱學習器結合形成強學習器來完成學習任務。Bagging 算法能夠有效減小預測的偏差和協(xié)方差，提高精度，避免過擬合［23］。本文使用ADRN 模型作為集成算法的基礎模型。

2.3 數據預處理

由于使用各類數據的空間分辨率不一致，存在MAIAC AOD 和ERA5 氣象數據的大區(qū)域覆蓋數據，也有類似HLH 和PM2.5濃度的點數據，為了便于制作模型訓練集，首先需對各類數據進行時空間匹配。遙感觀測數據中空間分辨率不足1 km 的參數，會使用雙線性插值法插值到1 km 分辨率；然后以MAIAC 觀測時間為基準，與其他氣象參數進行時間匹配；AQI 站點的PM2.5數據則與1 km 網格上的最近距離點進行匹配。匹配完成后進行標準化去除特征之間的相關性并使特征具有相似分布，提高模型訓練效率。

3 結果與分析

3.1 MAIAC AOD 填補

AOD 填補模型的輸入包括經度、緯度、地表高程、時間、地表濕度、溫度和壓強、東西風速、南北風速、邊界層高度、GEOS AOD 等特征，3 年時間的MAIAC 數據量太大，因此以周為單位進行數據集采樣和模型訓練。使用的深度學習框架為PYTORCH，訓練使用NVIDIA TITAN XP 顯卡，通過網格搜索尋找模型適合的超參數，最終ADRN模型參數設定如下：學習率為0.001，批量大小為512，訓練100 輪；采用學習率衰減策略，當監(jiān)測到當前損失函數在10 個輪次內沒有下降時，學習率縮小至1/2，以減小訓練時的參數震蕩，使得模型結果更接近最優(yōu)解；使用Adam 優(yōu)化器進行參數更新，在學習過程中自適應地迭代更新學習率和權重，相比隨機梯度下降有更好的效果。訓練完成后，模型的驗證結果R2為0.97。AOD 填補效果如圖4 所示，圖4（a）為實際MAIAC AOD 觀測結果，由于云層與地表因素影響，有較多區(qū)域無法反演得到AOD，存在大范圍空缺，圖4（b）、圖4（c）為模型填補后的AOD 分布。通過比較可以看出，填補的AOD 與原始MAIAC AOD 分布一致，且數據無明顯突變，分布連續(xù)，該填補模型能夠有效提高被動遙感數據的可用性。

圖4 2017 年1 月13 日MAIAC AOD 填補效果Fig.4 MAIAC AOD interpolation result on January 13th，2017

3.2 PM2.5濃度逐小時估算

在AOD 填補模型的基礎上，將填補后的MAIAC AOD 作為特征，與地基站點的PM2.5數據進行匹配后結合其他氣象數據形成PM2.5的訓練集，并利用機器學習算法與集成學習進行訓練測試。MLP 模型訓練超參數如下：3 層隱藏層，節(jié)點數分別為12、24 和12，迭代次數1 000 次，學習率為0.001，激活函數為RELU，批量大小為512，使用Adam 優(yōu)化器。SVR 模型訓練超參數如下：使用RBF 核函數，gamma 為0.06，訓練容差為0.001。RF模型訓練超參數為：最大隨機采樣率為50%，即每次隨機采樣樣本數占總訓練集的50%，最大樹深度20，決策樹數量500。集成學習的流程如下：首先對原始訓練數據集進行有放回地隨機采樣得到100 個樣本集，之后基于每個樣本集訓練基礎模型，得到100 個相互獨立的基學習器；其次進行實際PM2.5濃度估算時，將所有模型的輸出取平均作為集成模型的估計結果，ADRN 基礎模型的超參數與AOD 填補模型一致。

模型的測試結果如圖5 所示，MLP、SVR、RF和ADRN 集成模型分別對應圖5（a）～圖5（d），相關性分別達到0.43、0.55、0.83 和0.87，在相關性和均方根誤差（Root Mean Square Error，RMSE）上，ADRN 集成模型的表現最為優(yōu)異。同時比較了使用BLH 和HLH 時模型的驗證效果，結果見表1，使用HLH 特征對各模型的PM2.5的濃度估算結果都有明顯改善，如SVR 模型的R2從0.42 提高到0.55，同時偏差絕對值從0.62 減小到0.12，而ADRN 集成模型的相關性與偏差等指標在2 種情況下都最為優(yōu)異。

表1 使用 BLH 和 HLH 特征的PM2.5濃度估計結果Tab.1 Estimation results of PM2.5 concentration with BLH/HLH features

圖5 PM2.5濃度模型估計散點Fig.5 PM2.5 concentration model estimation

續(xù)圖5 PM2.5濃度模型估計散點Continuous Fig.5 PM2.5 concentration model estimation

訓練完成的ADRN 集成模型可用于估算研究區(qū)域在該3.5 a 時間內日間任意小時的PM2.5濃度分布，圖6（a）～圖6（d）展示了杭州市區(qū)2019 年12 月12 日在8：00、11：00、14：00 以及17：00 這幾個時間點的濃度分布估算結果，圖中的箭頭代表該時刻的風向。早上8：00 時，正值早高峰，交通排放增加，PM2.5高濃度地區(qū)集中在市中心并沿主要道路分布，最高濃度超過了100 μg/m3，此時整體風向為東北風；中午11：00 時，PM2.5濃度有一定程度的下降，且東部地區(qū)下降幅度較大；14：00 時，風向為偏東風，PM2.5平均濃度繼續(xù)下降，且東北部PM2.5濃度相比于西南角下降更多，東北風將高濃度地區(qū)的污染氣體攜帶到西部；17：00 時，城區(qū)的PM2.5濃度保持在50 μg/m3左右。PM2.5的濃度分布與日間交通狀況一致，且在交通干道周邊濃度較高，說明日間交通排放是城市顆粒物的主要來源之一。進一步分析PM2.5濃度的日間每小時及季節(jié)變化，如圖6（e）、圖6（f）所示。

圖6 杭州市區(qū)2019 年12 月12 日估算PM2.5濃度的空間分布及濃度變化時序Fig.6 Distribution of the estimated PM2.5 concentration in Hangzhou on December 12th，2019 and the concentration time-series

觀察到，日間不同時間的PM2.5濃度有明顯的變化，早上8～10 點濃度最高，之后逐漸下降，下午5 點后又有小幅上升，與早晚高峰時間匹配。同時觀察到，模型估計的PM2.5相較于站點實測數據更低，出現該現象的主要原因是AQI 監(jiān)測點基本分布在居民區(qū)及道路周邊，而居民區(qū)和道路周邊由于人類活動導致PM2.5濃度相對更高，而模型估算結果考慮的是區(qū)域內整體的PM2.5濃度，因此實測濃度比模型估算濃度高具有可解釋性。季節(jié)變化上，在冬季，12 月至來年2 月的濃度最高，夏季濃度最低。冬季空氣干燥，降水較少，且存在逆溫現象不利于污染氣體排放，同時冬季燃煤量和汽車尾氣排放加劇，這些因素共同導致冬季的PM2.5濃度居高不下；而夏季氣旋活動較頻繁，降水較多，有利于PM2.5的擴散和清除［24］。以上分析結果表明，以小時為精度的PM2.5濃度監(jiān)測對空氣質量的正確分級和有效管控有重要意義［25］。

4 結束語

本文提出了基于主被動結合和深度學習算法進行日間逐小時PM2.5濃度反演的方法：首先將從MPL 剖面計算得到的HLH 參數代替?zhèn)鹘y(tǒng)的BLH參數進行標準化；其次對MAIAC AOD、氣象參數與PM2.5濃度的相關性進行分析，改進了PM2.5逐時估算的準確性；最后基于MAIAC AOD 填補結果，建立了高時空分辨率的PM2.5濃度估算模型。通過多種模型的對比試驗，分析了HLH 對PM2.5濃度估算帶來的性能提升。以杭州城區(qū)為例，對該地區(qū)PM2.5濃度的時空變化進行了研究分析。

本研究對在監(jiān)測偏遠地區(qū)空氣質量與擴大PM2.5探測的時空覆蓋范圍有重要價值，而獲得高精度、高時空分辨率的空氣質量相關數據，是進行氣象、氣候、環(huán)境影響等諸多研究的必要基礎，也是對污染源進行準確定位、對排放量進行長期觀測，為有效的排查污染成因和制定治理政策提供參考的重要支撐［26］。在此研究基礎上，后續(xù)可結合長短期記憶神經網絡等循環(huán)神經網絡進行PM2.5濃度預測等工作。