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

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

0 引言

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

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

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

1.1 研究區(qū)域

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

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

1.2 地基數(shù)據(jù)

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

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

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

1.3 AOD 與氣象數(shù)據(jù)

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

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

2 研究方法與模型

2.1 機器學(xué)習(xí)模型

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

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

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

2.2 ADRN 與集成學(xué)習(xí)

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

圖3 ADRN 模型結(jié)構(gòu)及集成學(xué)習(xí)流程Fig.3 ADRN model structure and the integrated learning process

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

2.3 數(shù)據(jù)預(yù)處理

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

3 結(jié)果與分析

3.1 MAIAC AOD 填補

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

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

3.2 PM2.5濃度逐小時估算

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

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

表1 使用 BLH 和 HLH 特征的PM2.5濃度估計結(jié)果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

訓(xùn)練完成的ADRN 集成模型可用于估算研究區(qū)域在該3.5 a 時間內(nèi)日間任意小時的PM2.5濃度分布，圖6（a）～圖6（d）展示了杭州市區(qū)2019 年12 月12 日在8：00、11：00、14：00 以及17：00 這幾個時間點的濃度分布估算結(jié)果，圖中的箭頭代表該時刻的風(fēng)向。早上8：00 時，正值早高峰，交通排放增加，PM2.5高濃度地區(qū)集中在市中心并沿主要道路分布，最高濃度超過了100 μg/m3，此時整體風(fēng)向為東北風(fēng)；中午11：00 時，PM2.5濃度有一定程度的下降，且東部地區(qū)下降幅度較大；14：00 時，風(fēng)向為偏東風(fēng)，PM2.5平均濃度繼續(xù)下降，且東北部PM2.5濃度相比于西南角下降更多，東北風(fēng)將高濃度地區(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相較于站點實測數(shù)據(jù)更低，出現(xiàn)該現(xiàn)象的主要原因是AQI 監(jiān)測點基本分布在居民區(qū)及道路周邊，而居民區(qū)和道路周邊由于人類活動導(dǎo)致PM2.5濃度相對更高，而模型估算結(jié)果考慮的是區(qū)域內(nèi)整體的PM2.5濃度，因此實測濃度比模型估算濃度高具有可解釋性。季節(jié)變化上，在冬季，12 月至來年2 月的濃度最高，夏季濃度最低。冬季空氣干燥，降水較少，且存在逆溫現(xiàn)象不利于污染氣體排放，同時冬季燃煤量和汽車尾氣排放加劇，這些因素共同導(dǎo)致冬季的PM2.5濃度居高不下；而夏季氣旋活動較頻繁，降水較多，有利于PM2.5的擴散和清除［24］。以上分析結(jié)果表明，以小時為精度的PM2.5濃度監(jiān)測對空氣質(zhì)量的正確分級和有效管控有重要意義［25］。

4 結(jié)束語

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

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