董 焱,許 丹*,鮑艷松,許夢婕,顧英杰 (1.南京信息工程大學(xué),氣象災(zāi)害預(yù)報預(yù)警與評估協(xié)同創(chuàng)新中心,氣象環(huán)境衛(wèi)星工程與應(yīng)用聯(lián)合實驗室,中國氣象局氣溶膠與云降水重點開放實驗室,江蘇 南京 210044;2.南京信息工程大學(xué)大氣物理學(xué)院,江蘇 南京 210044)
大氣顆粒物濃度PM2.5是對人類健康影響較大的大氣污染指標(biāo),因為其顆粒直徑較小,粒子可以從肺部進入血管,同時有研究表明PM2.5微粒上附著的有害的物質(zhì)及重金屬可能進一步危害人類健康,目前PM2.5濃度的研究廣泛受到關(guān)注[1-2].氣溶膠光學(xué)厚度(AOD)是衛(wèi)星和地面遙感大氣氣溶膠容易獲得的光學(xué)參數(shù),通過衛(wèi)星還能得到區(qū)域氣溶膠分布特征.該參數(shù)表示整層大氣(從地面到大氣層頂)氣溶膠的削光量,是氣溶膠消光系數(shù)在垂直方向上的積分.由于氣溶膠本身對光路存在丁達(dá)爾效應(yīng),因此對大氣輻射和能量的收支平衡起著重要的作用,從而導(dǎo)致氣候變化的影響[3-4].因此大氣氣溶膠的定量是衛(wèi)星遙感的重要參數(shù)之一[5-6].
衛(wèi)星儀器探測的每一個像元都是反映整層大氣柱的氣溶膠光學(xué)特性,而大氣邊界層是人類活動的主要范圍.通過衛(wèi)星儀器(一般是可見光波段)估算近地面PM2.5濃度時,當(dāng)高層氣溶膠含量較高時,兩者的相關(guān)性會受到影響.粒子的吸濕作用如華東地區(qū)工業(yè)排放導(dǎo)致的粉煤灰以及未燃燒的碳在大氣中凝結(jié)水汽[7],導(dǎo)致AOD 受相對濕度影響較大.而顆粒物濃度結(jié)果不受相對濕度的影響.因此,AOD 與PM2.5在通常情況下存在非線性相關(guān)性[8].
利用衛(wèi)星傳感器在不同波段下觀測氣溶膠的光學(xué)特性,來估算觀測區(qū)域的PM2.5[9].其觀測結(jié)果具有高時效性、監(jiān)測成本低、且監(jiān)測范圍廣等優(yōu)勢[10-11].目前,關(guān)于此方面的工作,主要分為數(shù)學(xué)統(tǒng)計物理模型、人工智能,以及化學(xué)模式模擬等方法.其中數(shù)學(xué)統(tǒng)計的方法開展的最多,包括回歸模型[12]、多元回歸模型[13-14]、地理加權(quán)回歸[15]等,在此基礎(chǔ)上建立了氣溶膠細(xì)模態(tài)光學(xué)厚度與大氣顆粒物PM2.5質(zhì)量濃度之間的回歸關(guān)系[16];以及氣溶膠反射率與大氣顆粒物PM2.5質(zhì)量濃度的線性回歸關(guān)系等[17].目前對于人工智能/機器學(xué)習(xí)法開展的研究不多,基于神經(jīng)網(wǎng)絡(luò)方法對PM2.5質(zhì)量濃度進行估算[18].非線性的機器學(xué)習(xí)的方法也進行了很多不同算法之間的嘗試,但是和數(shù)學(xué)統(tǒng)計的方法一樣都需要較多的數(shù)據(jù)進行支撐,對數(shù)據(jù)的質(zhì)量的要求較高.而利用化學(xué)模式進行模擬,其依賴于模式本身.而物理機制建立模型的方法具有物理意義、通用型強、業(yè)務(wù)維護成本較低等特點.
由于我國國土面積遼闊、區(qū)域人口密度不協(xié)調(diào)、部分地區(qū)自然條件惡劣等一系列外在因素,想實現(xiàn)地面檢測網(wǎng)成為了一項具有挑戰(zhàn)性的任務(wù)[19].如今,華北地區(qū)形成了較為密集的觀測網(wǎng),但地面監(jiān)測網(wǎng)需要協(xié)調(diào)較多儀器,單個儀器之間的校準(zhǔn)定標(biāo)都存在一定差異.對于PM2.5的濃度局地趨勢延判具有較高難度.而通過衛(wèi)星觀測反演可以更好的解決上述問題.
利用FY4A/AGRI 儀器反演得到的AOD 數(shù)據(jù),以及對應(yīng)的地面相對濕度數(shù)據(jù),估算近地面大氣顆粒物PM2.5濃度.根據(jù)大氣氣溶膠消光的物理機制和垂直分布特征,運用大氣輻射傳輸模型6S(Second Simulation of the Satellite Signal in the Solar Spectrum)能見度經(jīng)驗轉(zhuǎn)換公式,計算得到氣溶膠標(biāo)高并對衛(wèi)星探測的整層大氣消光系數(shù)進行垂直訂正,估算近地面PM2.5質(zhì)量濃度數(shù).
FY4A 是中國發(fā)射的第二代靜止氣象衛(wèi)星,2016年12 月11 日發(fā)射升空,并在2017 年9 月25 日正式投入用戶使用[20].其搭載的掃描輻射成像儀躋身于世界靜止軌道成像儀最先進行列[21](Advanced Geostationary Radiation Imager,AGRI).AGRI 每15min 生成一副全圓盤影像觀測,且一共擁有14 個通道,其中兩個可見光通道(紅光),藍(lán)光以及近紅外、短波紅外、中波紅外和熱紅外通道等[22].除可見光通道外, AGRI 的空間分辨率在4km,這有利于多通道數(shù)據(jù)的計算,為利用數(shù)據(jù)進行同化和反演提供了較大的便利[23].
本文基于團隊自主研發(fā)的FY4A/AGRI 氣溶膠光學(xué)厚度數(shù)據(jù)集進行分析研究.該AOD 數(shù)據(jù)集采用暗像元法,利用AGRI 前6 個通道,基于6S 模型建立查算表反演得到,詳細(xì)原理和驗證過程另有文章闡明,在此不再贅述.本文利用此數(shù)據(jù)集,對PM2.5濃度進行估算.
搭載在NPP 衛(wèi)星上共有五個傳感器,其中包括VIIRS(Visible infrared Imaging Radiometer)即可見光紅外成像輻射儀.在2011 年10 月28 日啟動,星下點空間分辨率為400m[24],VIIRS 相對于MODIS 具有更高的分辨率[25].
圖1 華北東部地區(qū)環(huán)境監(jiān)測站點分布Fig.1 Map showing the location of the 180air quality monitoring stations
采用地面的相對濕度數(shù)據(jù)參與近地面顆粒物PM2.5的估算.同時,采用了地面顆粒物PM2.5數(shù)據(jù)來進行PM2.5估算的檢驗對比.利用美國METONE 公司所生產(chǎn)的手持式空氣塵埃粒子計數(shù)器進行近地面PM2.5測量.其數(shù)據(jù)覆蓋了京津冀和山東地區(qū).數(shù)據(jù)來源于各個省市的環(huán)境監(jiān)測中心180 個地面站點(圖1)信息,包括小時平均相對濕度數(shù)據(jù)和對應(yīng)的顆粒物數(shù)據(jù),如圖一所示.采用2018 年4 月~2019 年1 月不同季節(jié)高污染天氣前后5d 的數(shù)據(jù)樣本進行與地面數(shù)據(jù)的精度檢驗,并提取2018 年4 月18 日重污染天氣個例經(jīng)預(yù)處理后與地面172 個站點進行個例分析對比.
在原理上,由于衛(wèi)星所得到的AOD 數(shù)據(jù),每一個像元所對應(yīng)的是衛(wèi)星探測到的整層氣溶膠光學(xué)厚度即大氣柱上的總消光系數(shù),但是近地面的顆粒物PM2.5濃度所指的是對流層以下甚至是干粒子的質(zhì)量濃度,并且大氣氣溶膠的吸濕性對其消光特性具有相當(dāng)大的影響[26].所以AOD和PM2.5兩者的相關(guān)性受大氣分布和大氣狀態(tài)等一系列因素的影響[27].本文通過對AOD-PM2.5建模的分析,將影響較大因子放入模型中參與計算,具體計算流程如下圖2 所示.
圖2 實驗流程Fig.2 Experimental flow chart
2.2.1 垂直訂正 假設(shè)不同的大氣層各邊界之間相互平行,則AOD 是對整個大氣柱的氣溶膠消光系數(shù)積分,可以通過公式(1)進行表示:
式中:τa(λ)為AOD、Ka為氣溶膠消光系數(shù)在波長為λ 高度為z 時的值.可以通過估算的方式得到Ka(λ,z)的函數(shù)表達(dá)式[28],具體如下:
式中:Ka,o(λ)是表示地表氣溶膠的消光系數(shù)即z=0; HA為氣溶膠標(biāo)高,m;氣溶膠標(biāo)高可以看成大氣邊界層高度(Boundary Layer Height,簡稱BLH)[29].因此,式(1)、式(2)結(jié)合可以得到:
由式(3)中可以發(fā)現(xiàn),當(dāng)Ka,o(λ)、τa(λ)和HA如果有兩個是已知的,則可以計算到第三個.由于本文基于業(yè)務(wù)化系統(tǒng)進行考慮,對于大氣邊界層高度數(shù)據(jù)的獲取利用再分析資料的時效性都無法實現(xiàn)實時的業(yè)務(wù)化系統(tǒng)監(jiān)測,因此對于如何利用已知的AOD 數(shù)據(jù)來估算氣溶膠標(biāo)高變成了難題.通過在大氣輻射傳輸模型6S 中有AOD 與能見度的經(jīng)驗轉(zhuǎn)換模型:
式中:V 為近地面能見度m[30].增加在水平方向上的瑞利散射和臭氧吸收所帶來的影響[31-32],將0.02 的對比閾值更改到了0.0416,得到公式(5),具體如下:
根據(jù)式(5)估算出氣溶膠標(biāo)高HA.再將HA帶入式(4)中,得到垂直訂正后的消光系數(shù)Ka,o(λ).
2.2.2 濕度訂正 由于PM2.5指的是干粒子的質(zhì)量濃度,需要對AOD 數(shù)據(jù)進行濕度訂正[33].因此,吸濕性增長因子可以表示為[34]:
式中:RH 為地表相對濕度; f(RH)為吸濕性增長因子;g 為經(jīng)驗擬合系數(shù).華北東部地區(qū)的經(jīng)驗擬合系數(shù)g一般取值為0.38[35].
因此,濕度訂正后的消光系數(shù)Ka,o(λ)可以表示為[36]:
因此,將上述式(8)與式(5)結(jié)合,得到近地面干氣溶膠消光系數(shù)Ka,0,Dry(λ),公式如下:
2.2.3 氣溶膠類型參數(shù)化 根據(jù)粒子的米散射理論,大氣消光系數(shù)Ka定義為[37]:
式中:Qext(m, r ,λ )為大氣氣溶膠消光效率,m、r 分別表示單個粒子的復(fù)折射指數(shù)、粒子尺度,n(r)表示數(shù)密度譜分布.粒子尺度的取值范圍為(0,x/2)μm.
顆粒物PMx在大氣中的質(zhì)量濃度可以表示為:
式中:ρ為對應(yīng)尺度的氣溶膠質(zhì)量密度,單位為g/m3,北京附近地區(qū),ρ取1.5g/m3[38].Hansen等[39]研究發(fā)現(xiàn),消光效率Qext和有效半徑reff可以表示為:
表1 MODIS 在不同氣溶膠類型的消光特性表Table 1 Extinction characteristics of MODIS in different aerosol types
因此,將式(9)、(10)、(11)、(12)整理:
式中:Ka為濕度訂正和垂直訂正之后的大氣氣溶膠消光系數(shù),X 取2.5.在MODIS 的二級數(shù)據(jù)MOD04產(chǎn)品中,將不同氣溶膠類型(Aerosol Type 即Aerosol_Type_Land)進行了對消光效率Qext和有效半徑reff的分類,如表1 所示,選取其中大陸型的消光效率與有效半徑作為本次實驗的參數(shù).
將不同衛(wèi)星的AOD 數(shù)據(jù)進行預(yù)處理,包含對個別數(shù)據(jù)文件的邊角不整齊所進行的統(tǒng)一切割;對數(shù)據(jù)的缺省進行相應(yīng)的質(zhì)量控制;對數(shù)據(jù)在程序中占用空間太大所進行的數(shù)據(jù)格式轉(zhuǎn)化以及數(shù)據(jù)保留小數(shù)點后幾位的選擇等.
將每一個衛(wèi)星像素點匹配上與其像素點距離最近的地面信息數(shù)據(jù),如式(14)所示:
式中: Dmin表示為兩點的最小距離;lonsat表示衛(wèi)星某個像素點的經(jīng)度數(shù)值;latsat表示衛(wèi)星某個像素點的緯度數(shù)值.同樣,longro表示地面測站某個像素點的經(jīng)度數(shù)值,latgro表示地面測站某個像素點的緯度數(shù)值.
在時間匹配上,根據(jù)不同衛(wèi)星數(shù)據(jù)的過境時間信息尋找就近的相對濕度數(shù)據(jù)來進行估算以及對應(yīng)地面顆粒物濃度數(shù)據(jù)進行結(jié)果的檢驗.數(shù)據(jù)總結(jié)發(fā)現(xiàn):MODIS 與VIIRS 都以16d 為一個飛行周期,且每天飛過實驗地區(qū)為10:00~11:00.因為本文所使用的地面信息站點的數(shù)據(jù)為小時平均,與衛(wèi)星數(shù)據(jù)進行時間匹配時間窗小于±1h.
圖3 為FY-4A/AGRI PM2.5 區(qū)域分布圖,取自2019 年10 月28 日11:00 數(shù)據(jù)個例.在山東省的西南地區(qū)以及東部的半島沿海地區(qū)存在PM2.5濃度較大的區(qū)域.同時,北京和天津的PM2.5濃度同樣較高.
圖3 FY-4A/AGRI PM2.5 區(qū)域分布Fig.3 FY-4A/AGRI PM2.5 regional distribution map
圖4 是基于AGRI 計算出的PM2.5結(jié)果與地面站點匹配的數(shù)據(jù)折線圖,取自2018年4月18日11:00數(shù)據(jù).利用AGRI 得到的PM2.5與地面觀測值整體趨勢一致.在110號站點附件AGRI都存在高值,其原因為卷云的干擾或是云邊界像素,導(dǎo)致云識別沒有將此像素剔除,由于卷云反射率與地面不同,因此衛(wèi)星計算值偏高.
圖4 AGRI 反演結(jié)果與地面結(jié)果對比Fig.4 Comparison of AGRI inversion results and ground results
由于AOD 具有光學(xué)特性,其散射和吸收也會因不同季節(jié)太陽直射點的變化而受到影響.對顆粒物濃度夏季<秋季/春季<冬季,春季和秋季在PM2.5污染程度大致相同[40].本文隨機選取了夏季2019 年9 月18 日、秋季2018 年10 月13 日、春季2018 年4 月18 日和冬季2019 年1 月10 日進行對比,如圖5 所示.四個季節(jié)相比地面站點顆粒物PM2.5的高值區(qū)與低值區(qū)兩個衛(wèi)星都有較好的表現(xiàn),其次相比于地面顆粒物站點數(shù)值,衛(wèi)星數(shù)據(jù)估算的結(jié)果使其PM2.5的整體趨勢顯得更加直觀.在秋季(Ⅰ),無論是地面站點還是利用衛(wèi)星資料估算的顆粒物PM2.5都在天津地區(qū)存在一個高值區(qū),在河北和站東北部沿海地區(qū)具有較高的PM2.5質(zhì)量濃度.其質(zhì)量濃度數(shù)趨勢沿著海灣向四周遞減.利用衛(wèi)星估算的PM2.5結(jié)果與地面站點探測數(shù)據(jù)結(jié)果吻合.在春季(Ⅱ),地面站點PM2.5數(shù)據(jù)與衛(wèi)星估算PM2.5結(jié)果在北京地區(qū)附近都存在高值.兩個衛(wèi)星的結(jié)果具有較高的一致性,與地面站點探測數(shù)據(jù)相統(tǒng)一.在冬季(Ⅲ)MODIS 與VIIRS 對比差異相對較大.因為算法不同,因此導(dǎo)致舍去的像元不同.在MODIS 數(shù)據(jù)結(jié)果中,河北大部分地區(qū)包括北京和天津地區(qū)都存在缺省值,主要是這些地區(qū)沒有達(dá)到MODIS 官方暗像元的級別要求從而在反演中將其剔除.MODIS 數(shù)據(jù)刪去了大量的像素點,但是其保留下來的像素點與VIIRS 數(shù)據(jù)依然存在較高的一致性.相比于地面站點河北南部以及山東北部的顆粒物高值地帶,這兩顆衛(wèi)星都有較好的表現(xiàn).在夏季(Ⅳ),山東地區(qū)從沿海到內(nèi)陸PM2.5質(zhì)量濃度都有逐步變大的趨勢,兩顆衛(wèi)星結(jié)果都與地面站點大致相符,趨勢一致,但是VIIRS 不如MODIS 結(jié)果明顯.其原因主要是基于儀器標(biāo)定的差異以及VIIRS AOD 網(wǎng)格數(shù)據(jù)精度較低所導(dǎo)致.在地面站點圖中,北京地區(qū)具有一個PM2.5的低值區(qū).這兩個衛(wèi)星相對于地面結(jié)果同樣具有相同的趨勢.可以初步判斷,此方法可以用在多種衛(wèi)星估算顆粒物PM2.5質(zhì)量濃度上,受季節(jié)的變化影響較小,模型具有較高的穩(wěn)定性.
根據(jù)上述時空匹配的方法,對AGRI AOD 數(shù)據(jù)進行了多個個例結(jié)果的運行,將地面顆粒物PM2.5站點數(shù)據(jù)與AGRI、MODIS 和VIIRS PM2.5數(shù)據(jù)進行統(tǒng)計和評價.其中包括對散點對比、樣本個數(shù)(N)、相關(guān)系數(shù)(R)、均方根誤差(RMSE)、平均偏差(ME)以及相對誤差(MRE)等幾個方面的精度評價.在時間匹配上,采用的是地面站點PM2.5質(zhì)量濃度小時平均數(shù)據(jù);空間匹配上,地面PM2.5站點數(shù)據(jù)與衛(wèi)星像素點值的距離范圍小于0.01 經(jīng)緯度.若地面站點或衛(wèi)星像素點存在缺省則此匹配點不參與驗證.在質(zhì)量控制中,剔除了衛(wèi)星數(shù)據(jù)中的缺測值(-9999)以及地面站點的漏報值(-999);同時,若衛(wèi)星估算值和地面站點值差的絕對值大于其三倍的均方根誤差(RMSE)進行剔除,以保證其它外在因素的干擾.
將地面測站結(jié)果與衛(wèi)星數(shù)據(jù)結(jié)果進行上述預(yù)處理后,結(jié)果如圖6 和表2 所示.選取春季2018 年4 月14~18 日的10:00 數(shù)據(jù);夏季為2018 年7 月13~18 日10:00;秋季和冬季分別選取2018 年10月13~17 日10:00 和2019 年1 月10~14 日10:00數(shù)據(jù),分別對不同衛(wèi)星儀器的數(shù)據(jù)進行了精度檢驗分析.
圖5 不同季節(jié)的多顆衛(wèi)星反演結(jié)果與地面結(jié)果比較Fig.5 Comparison of satellite retrieval results and ground results in different seasons
圖6 四個季節(jié)不同AOD 數(shù)據(jù)估算近地面PM2.5 結(jié)果與地面測站PM2.5 結(jié)果散點對比Fig.6 Scatter comparison of PM2.5 results from different satellite instruments in four seasons and PM2.5 results from ground stations
如圖7 所示,采用中值濾波的方法對平均值進行平滑處理.其中,夏季(b)PM2.5平均濃度在京津冀地區(qū)達(dá)到了60μg/m3.其他3 個季節(jié)的平均濃度最大值約在80~100μg/m3之間.因此,夏季平均濃度相較于其他三個季節(jié)較低.相反,秋冬兩季,北京和天津地區(qū)的PM2.5平均濃度值為90~100μg/m3,比其他地區(qū)高,冬季更為明顯.并且除極值區(qū)外,冬季的低值區(qū)也有近50μg/m3,相較于其他3 個季節(jié)的低值區(qū)也相對偏高.進一步說明,夏季的顆粒物濃度普遍低,而冬季的顆粒物濃度普遍較高.
在以下4 個季節(jié)中AGRI 與地面數(shù)據(jù)相比都有相對較好的結(jié)果,相關(guān)系數(shù)R 都在0.8 以上,進一步反映出AGRI 的穩(wěn)定性.AGRI 在四個季節(jié)相對誤差(MRE)都小于20%,其中夏季達(dá)到17.24%.由表2 可知,AGRI 在冬季和春季的精度檢驗結(jié)果相對較好,夏季的精度普遍偏低[41].原因主要有兩點,其一是因為夏季的顆粒物濃度相對于其他三個季節(jié)濃度相對較低,衛(wèi)星儀器的靈敏程度對測量濃度較低的PM2.5偏差較大.其二是因為夏季的太陽高度角的變化使太陽直射點北移,導(dǎo)致暗像元算法受到季節(jié)的影響從而影響了暗像元識別的等級.但上述兩點仍需要進一步驗證.AGRI 的結(jié)果明顯優(yōu)于VIIRS 和MODIS 的結(jié)果.AGRI 在4 個季節(jié)中的相對誤差都小于VIIRS 和MODIS,AGRI 與地面的數(shù)據(jù)更加接近,AGRI 的穩(wěn)定性更好.在顆粒物濃度較大的春秋兩季,AGRI 的結(jié)果總體看也略優(yōu).夏季顆粒物濃度低于其他3 個季節(jié).因此夏季的散點在低值區(qū)普遍較多,在相關(guān)性上也略微偏低.但其結(jié)果仍不亞于MODIS 和VIIRS 的夏季精度檢驗結(jié)果.
圖7 AGRI PM2.5 季節(jié)區(qū)域分布結(jié)果Fig.7 Seasonal regional distribution of AGRI
表2 不同儀器在不同季節(jié)的精度驗證表Table 2 Accuracy verification table of different instruments in different seasons
本文采用FY-4A 的AOD 數(shù)據(jù)反演華北地區(qū)PM2.5濃度,基于輻射傳輸模型進行垂直訂正,同時結(jié)合RH數(shù)據(jù)進行濕度訂正.通過分析高污染天氣期間PM2.5濃度反演結(jié)果,并與VIIRS/MODIS 反演結(jié)果對比.
4.1 此方法區(qū)域個例的計算時間在3min 之內(nèi),提高了時間成本,為利用AGRI 實現(xiàn)顆粒物高濃度實時預(yù)警提供了可能.
4.2 AGRI 的結(jié)果在精度上不亞于 MODIS 和VIIRS,均方根誤差和相對偏差結(jié)果在4 個季節(jié)中都較優(yōu).進一步突出了AGRI 的穩(wěn)定性較高.
4.3 PM2.5觀測值與地面站點數(shù)值的相關(guān)性季節(jié)變化特征明顯.秋冬兩季在京津冀地區(qū)PM2.5平均濃度較高,而夏季整體濃度較低.在精度上,冬季相關(guān)系數(shù)略高于夏季,其影響因素仍需進一步驗證.