基于AGRI 顆粒物濃度遙感反演及季節變化分析

董 焱,許 丹*,鮑艷松,許夢婕,顧英杰 (1.南京信息工程大學,氣象災害預報預警與評估協同創新中心,氣象環境衛星工程與應用聯合實驗室,中國氣象局氣溶膠與云降水重點開放實驗室,江蘇 南京 210044；2.南京信息工程大學大氣物理學院,江蘇 南京 210044)

大氣顆粒物濃度PM2.5是對人類健康影響較大的大氣污染指標,因為其顆粒直徑較小,粒子可以從肺部進入血管,同時有研究表明PM2.5微粒上附著的有害的物質及重金屬可能進一步危害人類健康,目前PM2.5濃度的研究廣泛受到關注[1-2].氣溶膠光學厚度(AOD)是衛星和地面遙感大氣氣溶膠容易獲得的光學參數,通過衛星還能得到區域氣溶膠分布特征.該參數表示整層大氣(從地面到大氣層頂)氣溶膠的削光量,是氣溶膠消光系數在垂直方向上的積分.由于氣溶膠本身對光路存在丁達爾效應,因此對大氣輻射和能量的收支平衡起著重要的作用,從而導致氣候變化的影響[3-4].因此大氣氣溶膠的定量是衛星遙感的重要參數之一[5-6].

衛星儀器探測的每一個像元都是反映整層大氣柱的氣溶膠光學特性,而大氣邊界層是人類活動的主要范圍.通過衛星儀器(一般是可見光波段)估算近地面PM2.5濃度時,當高層氣溶膠含量較高時,兩者的相關性會受到影響.粒子的吸濕作用如華東地區工業排放導致的粉煤灰以及未燃燒的碳在大氣中凝結水汽[7],導致AOD 受相對濕度影響較大.而顆粒物濃度結果不受相對濕度的影響.因此,AOD 與PM2.5在通常情況下存在非線性相關性[8].

利用衛星傳感器在不同波段下觀測氣溶膠的光學特性,來估算觀測區域的PM2.5[9].其觀測結果具有高時效性、監測成本低、且監測范圍廣等優勢[10-11].目前,關于此方面的工作,主要分為數學統計物理模型、人工智能,以及化學模式模擬等方法.其中數學統計的方法開展的最多,包括回歸模型[12]、多元回歸模型[13-14]、地理加權回歸[15]等,在此基礎上建立了氣溶膠細模態光學厚度與大氣顆粒物PM2.5質量濃度之間的回歸關系[16];以及氣溶膠反射率與大氣顆粒物PM2.5質量濃度的線性回歸關系等[17].目前對于人工智能/機器學習法開展的研究不多,基于神經網絡方法對PM2.5質量濃度進行估算[18].非線性的機器學習的方法也進行了很多不同算法之間的嘗試,但是和數學統計的方法一樣都需要較多的數據進行支撐,對數據的質量的要求較高.而利用化學模式進行模擬,其依賴于模式本身.而物理機制建立模型的方法具有物理意義、通用型強、業務維護成本較低等特點.

由于我國國土面積遼闊、區域人口密度不協調、部分地區自然條件惡劣等一系列外在因素,想實現地面檢測網成為了一項具有挑戰性的任務[19].如今,華北地區形成了較為密集的觀測網,但地面監測網需要協調較多儀器,單個儀器之間的校準定標都存在一定差異.對于PM2.5的濃度局地趨勢延判具有較高難度.而通過衛星觀測反演可以更好的解決上述問題.

利用FY4A/AGRI 儀器反演得到的AOD 數據,以及對應的地面相對濕度數據,估算近地面大氣顆粒物PM2.5濃度.根據大氣氣溶膠消光的物理機制和垂直分布特征,運用大氣輻射傳輸模型6S(Second Simulation of the Satellite Signal in the Solar Spectrum)能見度經驗轉換公式,計算得到氣溶膠標高并對衛星探測的整層大氣消光系數進行垂直訂正,估算近地面PM2.5質量濃度數.

1 數據來源和研究區域

1.1 FY4A/AGRI AOD 數據

FY4A 是中國發射的第二代靜止氣象衛星,2016年12 月11 日發射升空,并在2017 年9 月25 日正式投入用戶使用[20].其搭載的掃描輻射成像儀躋身于世界靜止軌道成像儀最先進行列[21](Advanced Geostationary Radiation Imager,AGRI).AGRI 每15min 生成一副全圓盤影像觀測,且一共擁有14 個通道,其中兩個可見光通道(紅光),藍光以及近紅外、短波紅外、中波紅外和熱紅外通道等[22].除可見光通道外, AGRI 的空間分辨率在4km,這有利于多通道數據的計算,為利用數據進行同化和反演提供了較大的便利[23].

本文基于團隊自主研發的FY4A/AGRI 氣溶膠光學厚度數據集進行分析研究.該AOD 數據集采用暗像元法,利用AGRI 前6 個通道,基于6S 模型建立查算表反演得到,詳細原理和驗證過程另有文章闡明,在此不再贅述.本文利用此數據集,對PM2.5濃度進行估算.

1.2 NPP/VIIRS AOD 數據

搭載在NPP 衛星上共有五個傳感器,其中包括VIIRS(Visible infrared Imaging Radiometer)即可見光紅外成像輻射儀.在2011 年10 月28 日啟動,星下點空間分辨率為400m[24],VIIRS 相對于MODIS 具有更高的分辨率[25].

1.3 PM2.5 數據和相對濕度數據

圖1 華北東部地區環境監測站點分布Fig.1 Map showing the location of the 180air quality monitoring stations

采用地面的相對濕度數據參與近地面顆粒物PM2.5的估算.同時,采用了地面顆粒物PM2.5數據來進行PM2.5估算的檢驗對比.利用美國METONE 公司所生產的手持式空氣塵埃粒子計數器進行近地面PM2.5測量.其數據覆蓋了京津冀和山東地區.數據來源于各個省市的環境監測中心180 個地面站點(圖1)信息,包括小時平均相對濕度數據和對應的顆粒物數據,如圖一所示.采用2018 年4 月～2019 年1 月不同季節高污染天氣前后5d 的數據樣本進行與地面數據的精度檢驗,并提取2018 年4 月18 日重污染天氣個例經預處理后與地面172 個站點進行個例分析對比.

2 區域PM2.5 濃度反演原理

2.1 輻射傳輸模型

在原理上,由于衛星所得到的AOD 數據,每一個像元所對應的是衛星探測到的整層氣溶膠光學厚度即大氣柱上的總消光系數,但是近地面的顆粒物PM2.5濃度所指的是對流層以下甚至是干粒子的質量濃度,并且大氣氣溶膠的吸濕性對其消光特性具有相當大的影響[26].所以AOD和PM2.5兩者的相關性受大氣分布和大氣狀態等一系列因素的影響[27].本文通過對AOD-PM2.5建模的分析,將影響較大因子放入模型中參與計算,具體計算流程如下圖2 所示.

圖2 實驗流程Fig.2 Experimental flow chart

2.2 垂直訂正與濕度訂正

2.2.1 垂直訂正 假設不同的大氣層各邊界之間相互平行,則AOD 是對整個大氣柱的氣溶膠消光系數積分,可以通過公式(1)進行表示:

式中:τa(λ)為AOD、Ka為氣溶膠消光系數在波長為λ 高度為z 時的值.可以通過估算的方式得到Ka(λ,z)的函數表達式[28],具體如下:

式中:Ka,o(λ)是表示地表氣溶膠的消光系數即z=0; HA為氣溶膠標高,m;氣溶膠標高可以看成大氣邊界層高度(Boundary Layer Height,簡稱BLH)[29].因此,式(1)、式(2)結合可以得到:

由式(3)中可以發現,當Ka,o(λ)、τa(λ)和HA如果有兩個是已知的,則可以計算到第三個.由于本文基于業務化系統進行考慮,對于大氣邊界層高度數據的獲取利用再分析資料的時效性都無法實現實時的業務化系統監測,因此對于如何利用已知的AOD 數據來估算氣溶膠標高變成了難題.通過在大氣輻射傳輸模型6S 中有AOD 與能見度的經驗轉換模型:

式中:V 為近地面能見度m[30].增加在水平方向上的瑞利散射和臭氧吸收所帶來的影響[31-32],將0.02 的對比閾值更改到了0.0416,得到公式(5),具體如下:

根據式(5)估算出氣溶膠標高HA.再將HA帶入式(4)中,得到垂直訂正后的消光系數Ka,o(λ).

2.2.2 濕度訂正 由于PM2.5指的是干粒子的質量濃度,需要對AOD 數據進行濕度訂正[33].因此,吸濕性增長因子可以表示為[34]:

式中:RH 為地表相對濕度; f(RH)為吸濕性增長因子;g 為經驗擬合系數.華北東部地區的經驗擬合系數g一般取值為0.38[35].

因此,濕度訂正后的消光系數Ka,o(λ)可以表示為[36]:

因此,將上述式(8)與式(5)結合,得到近地面干氣溶膠消光系數Ka,0,Dry(λ),公式如下:

2.2.3 氣溶膠類型參數化 根據粒子的米散射理論,大氣消光系數Ka定義為[37]:

式中:Qext(m, r ,λ )為大氣氣溶膠消光效率,m、r 分別表示單個粒子的復折射指數、粒子尺度,n(r)表示數密度譜分布.粒子尺度的取值范圍為(0,x/2)μm.

顆粒物PMx在大氣中的質量濃度可以表示為:

式中:ρ為對應尺度的氣溶膠質量密度,單位為g/m3,北京附近地區,ρ取1.5g/m3[38].Hansen等[39]研究發現,消光效率Qext和有效半徑reff可以表示為:

表1 MODIS 在不同氣溶膠類型的消光特性表Table 1 Extinction characteristics of MODIS in different aerosol types

因此,將式(9)、(10)、(11)、(12)整理:

式中:Ka為濕度訂正和垂直訂正之后的大氣氣溶膠消光系數,X 取2.5.在MODIS 的二級數據MOD04產品中,將不同氣溶膠類型(Aerosol Type 即Aerosol_Type_Land)進行了對消光效率Qext和有效半徑reff的分類,如表1 所示,選取其中大陸型的消光效率與有效半徑作為本次實驗的參數.

2.3 星-地數據時空匹配

將不同衛星的AOD 數據進行預處理,包含對個別數據文件的邊角不整齊所進行的統一切割;對數據的缺省進行相應的質量控制;對數據在程序中占用空間太大所進行的數據格式轉化以及數據保留小數點后幾位的選擇等.

將每一個衛星像素點匹配上與其像素點距離最近的地面信息數據,如式(14)所示:

式中: Dmin表示為兩點的最小距離;lonsat表示衛星某個像素點的經度數值;latsat表示衛星某個像素點的緯度數值.同樣,longro表示地面測站某個像素點的經度數值,latgro表示地面測站某個像素點的緯度數值.

在時間匹配上,根據不同衛星數據的過境時間信息尋找就近的相對濕度數據來進行估算以及對應地面顆粒物濃度數據進行結果的檢驗.數據總結發現:MODIS 與VIIRS 都以16d 為一個飛行周期,且每天飛過實驗地區為10:00～11:00.因為本文所使用的地面信息站點的數據為小時平均,與衛星數據進行時間匹配時間窗小于±1h.

3 結果與討論

3.1 AGRI 個例結果展示

圖3 為FY-4A/AGRI PM2.5 區域分布圖,取自2019 年10 月28 日11:00 數據個例.在山東省的西南地區以及東部的半島沿海地區存在PM2.5濃度較大的區域.同時,北京和天津的PM2.5濃度同樣較高.

圖3 FY-4A/AGRI PM2.5 區域分布Fig.3 FY-4A/AGRI PM2.5 regional distribution map

圖4 是基于AGRI 計算出的PM2.5結果與地面站點匹配的數據折線圖,取自2018年4月18日11:00數據.利用AGRI 得到的PM2.5與地面觀測值整體趨勢一致.在110號站點附件AGRI都存在高值,其原因為卷云的干擾或是云邊界像素,導致云識別沒有將此像素剔除,由于卷云反射率與地面不同,因此衛星計算值偏高.

圖4 AGRI 反演結果與地面結果對比Fig.4 Comparison of AGRI inversion results and ground results

3.2 不同季節各衛星的對比

由于AOD 具有光學特性,其散射和吸收也會因不同季節太陽直射點的變化而受到影響.對顆粒物濃度夏季＜秋季/春季＜冬季,春季和秋季在PM2.5污染程度大致相同[40].本文隨機選取了夏季2019 年9 月18 日、秋季2018 年10 月13 日、春季2018 年4 月18 日和冬季2019 年1 月10 日進行對比,如圖5 所示.四個季節相比地面站點顆粒物PM2.5的高值區與低值區兩個衛星都有較好的表現,其次相比于地面顆粒物站點數值,衛星數據估算的結果使其PM2.5的整體趨勢顯得更加直觀.在秋季(Ⅰ),無論是地面站點還是利用衛星資料估算的顆粒物PM2.5都在天津地區存在一個高值區,在河北和站東北部沿海地區具有較高的PM2.5質量濃度.其質量濃度數趨勢沿著海灣向四周遞減.利用衛星估算的PM2.5結果與地面站點探測數據結果吻合.在春季(Ⅱ),地面站點PM2.5數據與衛星估算PM2.5結果在北京地區附近都存在高值.兩個衛星的結果具有較高的一致性,與地面站點探測數據相統一.在冬季(Ⅲ)MODIS 與VIIRS 對比差異相對較大.因為算法不同,因此導致舍去的像元不同.在MODIS 數據結果中,河北大部分地區包括北京和天津地區都存在缺省值,主要是這些地區沒有達到MODIS 官方暗像元的級別要求從而在反演中將其剔除.MODIS 數據刪去了大量的像素點,但是其保留下來的像素點與VIIRS 數據依然存在較高的一致性.相比于地面站點河北南部以及山東北部的顆粒物高值地帶,這兩顆衛星都有較好的表現.在夏季(Ⅳ),山東地區從沿海到內陸PM2.5質量濃度都有逐步變大的趨勢,兩顆衛星結果都與地面站點大致相符,趨勢一致,但是VIIRS 不如MODIS 結果明顯.其原因主要是基于儀器標定的差異以及VIIRS AOD 網格數據精度較低所導致.在地面站點圖中,北京地區具有一個PM2.5的低值區.這兩個衛星相對于地面結果同樣具有相同的趨勢.可以初步判斷,此方法可以用在多種衛星估算顆粒物PM2.5質量濃度上,受季節的變化影響較小,模型具有較高的穩定性.

3.3 反演結果精度檢驗

根據上述時空匹配的方法,對AGRI AOD 數據進行了多個個例結果的運行,將地面顆粒物PM2.5站點數據與AGRI、MODIS 和VIIRS PM2.5數據進行統計和評價.其中包括對散點對比、樣本個數(N)、相關系數(R)、均方根誤差(RMSE)、平均偏差(ME)以及相對誤差(MRE)等幾個方面的精度評價.在時間匹配上,采用的是地面站點PM2.5質量濃度小時平均數據;空間匹配上,地面PM2.5站點數據與衛星像素點值的距離范圍小于0.01 經緯度.若地面站點或衛星像素點存在缺省則此匹配點不參與驗證.在質量控制中,剔除了衛星數據中的缺測值(-9999)以及地面站點的漏報值(-999);同時,若衛星估算值和地面站點值差的絕對值大于其三倍的均方根誤差(RMSE)進行剔除,以保證其它外在因素的干擾.

將地面測站結果與衛星數據結果進行上述預處理后,結果如圖6 和表2 所示.選取春季2018 年4 月14～18 日的10:00 數據;夏季為2018 年7 月13～18 日10:00;秋季和冬季分別選取2018 年10月13～17 日10:00 和2019 年1 月10～14 日10:00數據,分別對不同衛星儀器的數據進行了精度檢驗分析.

圖5 不同季節的多顆衛星反演結果與地面結果比較Fig.5 Comparison of satellite retrieval results and ground results in different seasons

圖6 四個季節不同AOD 數據估算近地面PM2.5 結果與地面測站PM2.5 結果散點對比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平均濃度在京津冀地區達到了60μg/m3.其他3 個季節的平均濃度最大值約在80～100μg/m3之間.因此,夏季平均濃度相較于其他三個季節較低.相反,秋冬兩季,北京和天津地區的PM2.5平均濃度值為90～100μg/m3,比其他地區高,冬季更為明顯.并且除極值區外,冬季的低值區也有近50μg/m3,相較于其他3 個季節的低值區也相對偏高.進一步說明,夏季的顆粒物濃度普遍低,而冬季的顆粒物濃度普遍較高.

在以下4 個季節中AGRI 與地面數據相比都有相對較好的結果,相關系數R 都在0.8 以上,進一步反映出AGRI 的穩定性.AGRI 在四個季節相對誤差(MRE)都小于20%,其中夏季達到17.24%.由表2 可知,AGRI 在冬季和春季的精度檢驗結果相對較好,夏季的精度普遍偏低[41].原因主要有兩點,其一是因為夏季的顆粒物濃度相對于其他三個季節濃度相對較低,衛星儀器的靈敏程度對測量濃度較低的PM2.5偏差較大.其二是因為夏季的太陽高度角的變化使太陽直射點北移,導致暗像元算法受到季節的影響從而影響了暗像元識別的等級.但上述兩點仍需要進一步驗證.AGRI 的結果明顯優于VIIRS 和MODIS 的結果.AGRI 在4 個季節中的相對誤差都小于VIIRS 和MODIS,AGRI 與地面的數據更加接近,AGRI 的穩定性更好.在顆粒物濃度較大的春秋兩季,AGRI 的結果總體看也略優.夏季顆粒物濃度低于其他3 個季節.因此夏季的散點在低值區普遍較多,在相關性上也略微偏低.但其結果仍不亞于MODIS 和VIIRS 的夏季精度檢驗結果.

圖7 AGRI PM2.5 季節區域分布結果Fig.7 Seasonal regional distribution of AGRI

表2 不同儀器在不同季節的精度驗證表Table 2 Accuracy verification table of different instruments in different seasons

4 結論

本文采用FY-4A 的AOD 數據反演華北地區PM2.5濃度,基于輻射傳輸模型進行垂直訂正,同時結合RH數據進行濕度訂正.通過分析高污染天氣期間PM2.5濃度反演結果,并與VIIRS/MODIS 反演結果對比.

4.1 此方法區域個例的計算時間在3min 之內,提高了時間成本,為利用AGRI 實現顆粒物高濃度實時預警提供了可能.

4.2 AGRI 的結果在精度上不亞于 MODIS 和VIIRS,均方根誤差和相對偏差結果在4 個季節中都較優.進一步突出了AGRI 的穩定性較高.

4.3 PM2.5觀測值與地面站點數值的相關性季節變化特征明顯.秋冬兩季在京津冀地區PM2.5平均濃度較高,而夏季整體濃度較低.在精度上,冬季相關系數略高于夏季,其影響因素仍需進一步驗證.