李莉，李長軍，趙海軍，葛瑞婷

(1.山東省氣象防災(zāi)減災(zāi)重點實驗室，山東 濟(jì)南 250031；2.臨沂市氣象局，山東 臨沂 276004；3.山東省氣象信息中心，山東 濟(jì)南 250031；4.淄博市氣象局，山東 淄博 255022)

引言

近年來由于經(jīng)濟(jì)高速發(fā)展，城市化進(jìn)程不斷推進(jìn)，環(huán)境污染問題日益嚴(yán)重，大氣污染事件頻發(fā)，引起了政府和社會各界的廣泛關(guān)注，大氣環(huán)境污染的監(jiān)測和治理成為大家關(guān)注的重點[1-2]。氣溶膠是影響大氣環(huán)境質(zhì)量的重要因素，在大氣中以固態(tài)或液態(tài)形式存在，通常指的是懸浮性顆粒。氣溶膠光學(xué)厚度 (aerosol optical depth，AOD)通常用來描述氣溶膠對光的削弱能力，指的是介質(zhì)的消光系數(shù)在垂直方向上的積分，可以計算大氣中的氣溶膠含量，并且能在一定程度上反映氣溶膠濃度以及空氣的污濁程度和污染情況[3-4]。評估氣溶膠光學(xué)厚度的值能用于反映氣溶膠粒子占比情況以及空氣是否清潔，在研究氣候狀況和環(huán)境污染情況方面有著重要的意義，開展氣溶膠光學(xué)厚度的地基觀測與遙感反演已成為大氣環(huán)境等領(lǐng)域的重要研究方向[5-7]。

傳統(tǒng)的地基觀測雖然數(shù)據(jù)精度高，但只能提供站點附近的局部氣溶膠信息，并且受人力物力的限制，目前的監(jiān)測站點稀疏，無法提供大的區(qū)域范圍的氣溶膠空間分布的細(xì)節(jié)信息，而且維護(hù)成本較高[8]。相比于地基監(jiān)測，衛(wèi)星遙感技術(shù)監(jiān)測范圍廣，空間覆蓋好，時間序列長，時效性強(qiáng)，成本低，是對地基觀測的有效補(bǔ)充[9]。近年來靜止衛(wèi)星由于空間分辨率顯著提高，被廣泛用于大氣污染物監(jiān)測，它的高度比較高，因而視野寬闊，可以對同一區(qū)域?qū)崿F(xiàn)連續(xù)觀測，和極軌衛(wèi)星相比，具有比較明顯的優(yōu)勢[10]。

國內(nèi)外已有不少關(guān)于氣溶膠的遙感反演研究。利用暗目標(biāo)像元法融合多種衛(wèi)星數(shù)據(jù)反演氣溶膠光學(xué)厚度，比如MODIS L1B級數(shù)據(jù)、高分四號衛(wèi)星數(shù)據(jù)[11-12]；并在此研究基礎(chǔ)上，又對氣溶膠光學(xué)厚度產(chǎn)品進(jìn)行更深一層的研究。關(guān)于利用靜止衛(wèi)星數(shù)據(jù)反演AOD方面：牛曉君等[13]提出一種新算法，即利用葵花衛(wèi)星可見光波段數(shù)據(jù)及地表反射率比值反演；關(guān)雷等[14]利用暗像元法基于高分四號衛(wèi)星數(shù)據(jù)，針對哈爾濱地區(qū)進(jìn)行了AOD反演。最近幾年也有針對我國風(fēng)云衛(wèi)星系列的研究成果，比如利用FY-2C靜止衛(wèi)星可見光資料[15]、FY-3A/MERSI(Medium Resolution Spectral Imager)數(shù)據(jù)[16]、FY-3C/MERSI數(shù)據(jù)[17]研究AOD的反演算法。綜上所述，目前各種衛(wèi)星傳感器反演AOD的經(jīng)驗豐富，但由于反演需要構(gòu)建本地化的氣溶膠模型，導(dǎo)致難以大范圍開展研究，所以針對FY-4A/AGRI(Advanced Geosynchronous Radiation Imager)傳感器數(shù)據(jù)反演AOD的研究仍然比較缺乏。鑒于此，本文以山東地區(qū)為研究區(qū)域，自定義氣溶膠模型，選取2021年兩次典型天氣污染過程利用FY-4A/AGRI數(shù)據(jù)進(jìn)行AOD反演，并將反演結(jié)果與鄰近時刻的氣溶膠光學(xué)厚度產(chǎn)品進(jìn)行細(xì)節(jié)表達(dá)、空間一致性的對比驗證以及定量分析，旨在探討適用于山東地區(qū)的自定義氣溶膠模型，以及新型國產(chǎn)衛(wèi)星FY-4A數(shù)據(jù)用于氣溶膠反演的可行性。

1 數(shù)據(jù)

1.1 FY-4A/AGRI衛(wèi)星遙感數(shù)據(jù)

FY-4A于2016年12月11日00：11(北京時，下同)發(fā)射升空，并于2018年初完成在軌測試，同年5月1日正式投入業(yè)務(wù)運行，是第二代靜止氣象衛(wèi)星的首發(fā)星。它位于104.7°E的赤道上空，采用三軸穩(wěn)定姿態(tài)控制，與風(fēng)云二號氣象衛(wèi)星的自旋穩(wěn)定方式相比，可顯著提高對地觀測效率。FY-4A/AGRI觀測數(shù)據(jù)可在風(fēng)云衛(wèi)星遙感數(shù)據(jù)服務(wù)網(wǎng)(http://satellite.nsmc.org.cn/PortalSite/Data/Satellite.aspx#)下載，下載格式包括中國區(qū)域和全圓盤兩種，分辨率有500 m、1 km、2 km、4 km[18-20]。研究利用FY-4A成像儀全圓盤2 km L1數(shù)據(jù)、L2級云產(chǎn)品和4 km 地球靜止軌道(geostationary orbit，GEO)數(shù)據(jù)，利用可見光的紅藍(lán)波段以及中紅外波段開展山東地區(qū)氣溶膠的反演研究工作。

1.2 MODIS地表反射率數(shù)據(jù)

地表反射率數(shù)據(jù)庫是氣溶膠反演的基礎(chǔ)，分辨率為500 m的8 d合成產(chǎn)品MOD09A1，選取8 d中質(zhì)量最好的數(shù)據(jù)合成得到地表反射率數(shù)據(jù)，可以有效消除短期內(nèi)的變化，并且考慮到云、氣溶膠等的影響及吸收、散射作用，最大程度上減少云覆蓋對地表反射率的影響[21]。MODIS數(shù)據(jù)的獲取渠道很多，很多國家都有自己的數(shù)據(jù)接收站點，本研究使用的MOD09A1(MODIS Terra Surface Reflectance 8-Day L3 Global 500 m SIN Grid)地表反射率產(chǎn)品，由美國航空航天局(National Aeronautics and Space Administration，NASA)Reverb數(shù)據(jù)中心(http://reverb.echo.nasa.gov)提供。

1.3 MODIS氣溶膠產(chǎn)品MOD04

MODIS的氣溶膠產(chǎn)品已經(jīng)比較成熟，穩(wěn)定運行20 a左右，并且已經(jīng)在全球范圍內(nèi)得到驗證，可以作為比較準(zhǔn)確的標(biāo)準(zhǔn)數(shù)據(jù)來驗證氣溶膠的反演結(jié)果。NASA官方發(fā)布的L2級近實時的氣溶膠產(chǎn)品MOD04_3K(MODIS Terra Aerosol 5-min L2 Swath 3 km)是本次用于驗證反演結(jié)果準(zhǔn)確性的數(shù)據(jù)之一，采用暗像元法進(jìn)行反演，它的空間分辨率為3 km、0.1°，能夠方便及時地獲得全球的海洋表面、陸地的大氣氣溶膠光學(xué)特性及其參數(shù)[22]。以上所有數(shù)據(jù)可在LAADS(Level 1 and Atmosphere Archive and Distribution System，http://ladsweb.nascom.nasa.gov/data/search.-html)搜索并下載。

1.4 Himawari-8氣溶膠光學(xué)厚度遙感產(chǎn)品

日本氣象廳的新一代靜止氣象衛(wèi)星Himawari-8于2014年10月成功發(fā)射，2015年7月7日成功通過測試開始投入正常使用，該衛(wèi)星具有觀測范圍廣、觀測周期短、空間分辨率高等優(yōu)點，它搭載的傳感器AHI(Advanced Himawari Imager)可實現(xiàn)10 min/次的高時間分辨率對地觀測。Himawari-8衛(wèi)星在140°E的上空，能夠觀測120°×120°的范圍，即80°E～20°W，60°N～60°S，空間分辨率為5 km×5 km，可以較好地觀測中國的大部分區(qū)域。本研究使用的是日本宇宙航空研究開發(fā)機(jī)構(gòu)(Japan Aerospace Exploration Agency，簡稱JAXA)官方發(fā)布的L2級氣溶膠光學(xué)厚度產(chǎn)品，該產(chǎn)品是YOSHIDA et al.[23]基于前人的研究基礎(chǔ)，在目標(biāo)函數(shù)中引入每個通道的權(quán)重，選擇用于氣溶膠反演的最佳通道，然后在陸地和海洋上設(shè)置候選氣溶膠模型，最后在 300～2 500 nm波長范圍內(nèi)每1 nm編制一次查找表，并使用每個傳感器的響應(yīng)函數(shù)對輻射進(jìn)行加權(quán)反演得到的。

1.5 其他反演所需數(shù)據(jù)

通過氣溶膠自動觀測網(wǎng)(Aerosol Robotic Network，AERONET)可以得到全球各個測站的氣溶膠光學(xué)厚度，再通過光學(xué)厚度以及其他資料的反演可以得到氣溶膠粒徑譜分布、單次散射反照率、吸收性光學(xué)厚度以及體積濃度等產(chǎn)品[24]。AERONET的氣溶膠產(chǎn)品可用于驗證與氣溶膠有關(guān)的衛(wèi)星產(chǎn)品以及實現(xiàn)與其他儀器的協(xié)同觀測等。為實現(xiàn)該目標(biāo)，AERONET 對儀器、測量、校準(zhǔn)和處理有標(biāo)準(zhǔn)化要求，使其在科學(xué)界有了很大的擴(kuò)展和廣泛的使用。本研究利用山東臨近站點徐州2013—2019年的L1.5級數(shù)據(jù)的粒子譜分布、復(fù)折射指數(shù)數(shù)據(jù)，計算對應(yīng)參數(shù)的季均值結(jié)果，把這些結(jié)果參數(shù)導(dǎo)入6S模式，從而建立自定義氣溶膠模型的反演查找表。

2 研究方法

2.1 氣溶膠光學(xué)厚度反演原理

在大氣遙感中，需要通過大氣輻射傳輸方程來求解氣溶膠光學(xué)厚度。假設(shè)地球表面為均勻的朗伯體，并且大氣在水平方向上的變化比較均勻，那么可見光通道衛(wèi)星傳感器接收到的相關(guān)輻射信息包括大氣散射、地表反射與大氣耦合兩個部分，衛(wèi)星接收到的大氣頂層表觀反射率可以表示為[25]：

(1)

式中：ρTOA(θs，θv，φ)是表觀反射率，又稱大氣頂部(top of atmosphere,TOA)反射率；ρ0(θs，θv，φ)是大氣程輻射，即大氣路徑輻射項等效反射率；S為大氣下界半球反射率；T(θs)、T(θv)是對應(yīng)方向上的大氣透過率，可合并寫為一項T；ρs(θs，θv，φ)是地表反射率，φ為相對方位角，θv為觀測天頂角，θs為太陽天頂角。其中，ρ0、T、S均為氣溶膠光學(xué)厚度的函數(shù)。衛(wèi)星遙感氣溶膠就是把這兩部分信息都估算出來，重構(gòu)TOA反射率的過程[26-27]。當(dāng)重構(gòu)出的TOA反射率最接近于衛(wèi)星觀測的TOA反射率時，即可得到最接近實際的ρ0、T、S，又由于這3個參數(shù)與AOD相關(guān)，得到這3個參數(shù)實際上就能得到AOD。所以，要反演出AOD，必須估算出地表反射率和3個大氣輻射參數(shù)，地表反射率主要通過暗目標(biāo)(Dark Target，DT)算法實現(xiàn)，3個大氣輻射參數(shù)的獲得主要通過6S輻射傳輸模型計算，構(gòu)建出AOD反演查找表。另外，由于有云條件下很難得到AOD，一般對有云像元不進(jìn)行反演，所以還需要對地表反射率產(chǎn)品數(shù)據(jù)進(jìn)行云檢測[28]，本文采用閾值法來進(jìn)行云檢測，利用下墊面和云在可見光、紅外波段的反射率差異進(jìn)行判斷，即當(dāng)0.55 μm波段的表觀反射率ρ0.55方差大于0.002 5時，則判定為云；當(dāng)1.61 μm波段表觀反射率ρ1.61大于0.25，判斷該像元為云；當(dāng)1.37 μm波段表觀反射率ρ1.37大于0.3時，則判斷該像元為云。具體反演流程如圖1所示。

圖1 氣溶膠光學(xué)厚度反演流程圖Fig.1 Flow chart of retrieving aerosol optical depth

綜上所述，利用遙感衛(wèi)星數(shù)據(jù)進(jìn)行氣溶膠光學(xué)厚度反演，主要有2個問題：一是確定符合當(dāng)?shù)貙嶋H情況的氣溶膠模型，另一個是剔除地表反射率的貢獻(xiàn)[29]。在進(jìn)行氣溶膠光學(xué)厚度反演的過程中，一般要先對反演區(qū)域的氣溶膠模型進(jìn)行假設(shè)，在此基礎(chǔ)上分別確定3個大氣輻射參數(shù)與AOD的對應(yīng)關(guān)系，再確定反演影像上的像元地表反射率。通過像元的地表反射率來確認(rèn)該像元的AOD值，從而實現(xiàn)反演。

2.2 地表反射率估算

反演過程中，在構(gòu)建地表反射率庫時，選取2018—2020年這3 a同一時期的地表反射率產(chǎn)品MOD09，該產(chǎn)品8 d合成。由于MODIS的地表反射率產(chǎn)品受云、云陰影、氣象要素、天氣條件的影響，以及每8 d才有一景影像，不能和FY-4A衛(wèi)星數(shù)據(jù)的時間一一對應(yīng)，如果直接用于AOD的反演會增加誤差。為了減小影響，假設(shè)大部分地物在短時間內(nèi)的地表反射率變化較小，由于MOD09是8 d合成，所以構(gòu)建數(shù)據(jù)集時將周期設(shè)定為3個月[30]。同時，為了避免受云影響嚴(yán)重，設(shè)置云覆蓋標(biāo)準(zhǔn)小于35%，如果某個時間段內(nèi)符合條件的影像較少或沒有，可以適當(dāng)放寬云覆蓋標(biāo)準(zhǔn)或選用臨近年份對應(yīng)時期的數(shù)據(jù)代替，從而保證每個周期內(nèi)都有足夠的影像參與數(shù)據(jù)集的構(gòu)建。對MODIS地表反射率產(chǎn)品進(jìn)行投影轉(zhuǎn)換、拼接、裁剪、校正、重采樣等預(yù)處理，用以獲取研究區(qū)規(guī)定時間內(nèi)的地表反射率信息。然后根據(jù)構(gòu)建的反射率庫確定其與FY-4A的表觀反射率信息之間的線性關(guān)系，將線性關(guān)系應(yīng)用于滿足暗像元的清潔像元，獲取地氣解耦的初步結(jié)果，結(jié)合觀測幾何等，進(jìn)一步判定大氣貢獻(xiàn)部分。根據(jù)自定義的氣溶膠模型以及衛(wèi)星拍攝時刻下的觀測幾何，讀入各個像元的觀測幾何，實現(xiàn)對不同像元在自定義氣溶膠模型下估算大氣部分，以此將地表反射率從衛(wèi)星探測到的表觀反射率信號中剝離，實現(xiàn)地氣分離。

2.3 自定義氣溶膠模型

自定義氣溶膠模型，也就是氣溶膠模型的本地化確定，是氣溶膠光學(xué)厚度反演的另一關(guān)鍵點。氣溶膠的基本組成包括沙塵性粒子、水溶性粒子、海洋性粒子和煤煙性粒子，不同粒子的大氣散射效應(yīng)不同，對大氣輻射傳輸?shù)挠绊懸膊煌?，而不同的氣溶膠模型中所包含的粒子種類及比例不同，物理特性和光學(xué)特性也就隨之不同，復(fù)折射指數(shù)、對太陽輻射的吸收和散射等也就不一樣，所以選用不同的氣溶膠類型，將導(dǎo)致反演的結(jié)果有較大的差異[31-32]。氣溶膠模型實際上就代表了不同類型的粒子，代表氣溶膠不同的組分及其所占的不同比率。山東地區(qū)地形復(fù)雜，既有山地丘陵，又有沿海城市，在春季受到沙塵影響，在夏秋季受到秸稈焚燒影響，在冬季則會受到燃煤取暖影響[33]，僅用一種氣溶膠模型不能滿足實際反演需要?；谏鲜隹紤]，結(jié)合山東地區(qū)實際情況選擇從AERONET網(wǎng)站下載山東臨近站點徐州2013—2019年的L1.5級數(shù)據(jù)，對山東區(qū)域進(jìn)行氣溶膠模型的本地化構(gòu)建，同時利用衛(wèi)星、地基觀測的表觀反射率數(shù)據(jù)和地表反射率數(shù)據(jù)，結(jié)合輻射傳輸方程，使用反復(fù)隨機(jī)抽樣方法來模擬氣溶膠模型。按季節(jié)對粗細(xì)模態(tài)半徑、粗細(xì)模態(tài)體積濃度、復(fù)折射指數(shù)進(jìn)行統(tǒng)計，得到對應(yīng)參數(shù)的均值及標(biāo)準(zhǔn)差范圍。在各參數(shù)的標(biāo)準(zhǔn)差范圍內(nèi)，以步長0.001結(jié)合輻射傳輸模型進(jìn)行迭代。在迭代時，將當(dāng)前迭代下的氣溶膠參數(shù)結(jié)果、地表反射率結(jié)果、氣溶膠光學(xué)厚度觀測值輸入輻射傳輸模型，計算出當(dāng)前季節(jié)下所有的理論表觀反射率和衛(wèi)星實際表觀反射率的殘差。所有迭代完成后，尋找最小殘差和對應(yīng)的氣溶膠模型迭代參數(shù)，將該參數(shù)作為對應(yīng)季節(jié)的氣溶膠模型參數(shù)。最終得到山東省季節(jié)性氣溶膠類型參數(shù)，如表1所示。

表1 山東省季節(jié)性氣溶膠類型參數(shù)

2.4 查找表的構(gòu)建

由于大氣輻射傳輸模型的計算速度比較緩慢，在反演前需要建立查找表。本研究選擇6S大氣輻射傳輸模型，在生成多維查找表的過程中，根據(jù)研究區(qū)域的地理特征、大氣條件，并且考慮山東地區(qū)地形、氣溶膠特性等因素，設(shè)置了適用于山東區(qū)域的參數(shù)，以此生成滿足反演區(qū)域的自定義氣溶膠類型的兩套紅藍(lán)波段查找表。根據(jù)完成地氣解耦后獲取的地表反射率以及FY-4A數(shù)據(jù)上取得的表觀反射率建立查找關(guān)系，通過交互式數(shù)據(jù)語言(interactive data language, IDL)將地氣解耦后獲取的地表反射率文件讀入，然后遍歷能夠反演氣溶膠光學(xué)厚度的像元，讀取其對應(yīng)的地表反射率以及觀測幾何等參數(shù)，以此為根據(jù)從查找表中獲取對應(yīng)臨近的氣溶膠光學(xué)厚度的值，再經(jīng)過線性插值后獲得0.47、0.66 μm兩個波段的氣溶膠光學(xué)厚度，對兩者求平均獲取0.55 μm的氣溶膠光學(xué)厚度，完成對研究區(qū)域逐像元的氣溶膠光學(xué)厚度計算。

大氣模型中的主要參數(shù)是臭氧和水汽含量，這兩個參數(shù)對用于氣溶膠反演的波段影響較小，在研究中不重點考慮，因此將6S模型中的大氣模型設(shè)置為中緯度冬季和中緯度夏季兩種模型，衛(wèi)星天頂角、太陽天頂角等其他輸入?yún)?shù)設(shè)置詳見表2。

表2 查找表具體參數(shù)設(shè)置

3 AOD反演結(jié)果與精度驗證

根據(jù)上述研究方法，準(zhǔn)備反演所需的地表反射率庫、查找表等數(shù)據(jù)，在此基礎(chǔ)上，選取山東地區(qū)2021年1月22—24日、10月28—30日兩次典型污染案例的AGRI數(shù)據(jù)進(jìn)行反演試驗，并將結(jié)果分別與MODIS、Himawari-8/AHI的氣溶膠光學(xué)厚度產(chǎn)品進(jìn)行對比分析。不同時刻的AGRI反演的氣溶膠結(jié)果，分別對應(yīng)其相同或臨近時刻的MODIS、Himawari-8/AHI氣溶膠光學(xué)厚度產(chǎn)品數(shù)據(jù)。目前使用這兩種產(chǎn)品進(jìn)行精度驗證分析的研究比較多，這兩種產(chǎn)品的算法比較權(quán)威，產(chǎn)品也比較成熟，應(yīng)用比較普遍，可以用來評價不同傳感器的氣溶膠反演結(jié)果。

3.1 FY-4A衛(wèi)星2021年1月22—24日反演結(jié)果

輸入6S查找表、預(yù)處理過的FY-4A全圓盤2 km L1、4 km GEO數(shù)據(jù)等，進(jìn)行反演，得到AOD結(jié)果。圖2a、b、c分別為1月22日12：00、23日12：00、24日12：00 AGRI的反演結(jié)果，圖中從紅到藍(lán)表示氣溶膠光學(xué)厚度逐漸減小，反演結(jié)果的空白區(qū)域是有云覆蓋或者地表反射率過高不滿足反演條件。

通過查閱氣象信息及相關(guān)資料可知，22—24日受不利氣象條件影響，山東省內(nèi)陸地區(qū)出現(xiàn)一次污染過程。22日，地面處于弱氣壓場，魯中、魯南地區(qū)受氣流輻合影響，擴(kuò)散條件較差，部分地區(qū)出現(xiàn)中至重度污染，其他地區(qū)受偏東、偏北氣流控制，擴(kuò)散條件稍好，魯西北部分地區(qū)出現(xiàn)輕至中度污染。23日，沿海地區(qū)受持續(xù)偏東氣流控制，空氣質(zhì)量較好，而魯西南地區(qū)有氣流輻合，濕度較大，擴(kuò)散條件較差，出現(xiàn)輕至重度污染。24日，受冷空氣和降雪共同影響，污染過程結(jié)束，全省空氣質(zhì)量恢復(fù)到良好水平。

圖2 2021年1月22—24日FY-4A反演AOD(色階)結(jié)果分布(a. 22日， b. 23日， c. 24日)Fig.2 Distribution of FY-4A retrieving AOD results (color scale) on 22 (a), 23 (b), and 24 (c) January 2021

通過圖2可以看出， 1月23日在魯西南地區(qū)監(jiān)測到輕至重度污染，在西南部存在分布較為密集的AOD高值區(qū)域，AOD值普遍在0.6以上。對比圖2a、b、c發(fā)現(xiàn)，魯西南部分地區(qū)22日出現(xiàn)中至重度污染，23日污染范圍擴(kuò)大，24日空氣質(zhì)量好轉(zhuǎn)，在圖上表示為顏色由黃變紅再變藍(lán)。圖2a和圖2c中，全省大部分區(qū)域AOD值較低，在圖上呈現(xiàn)深藍(lán)色，AOD高值區(qū)域零星分布在魯西南區(qū)域(圖2a)，山東省大部分區(qū)域AOD值較小(圖2c)。據(jù)此可以推測，此次空氣污染過程已經(jīng)結(jié)束。綜上所述，基于AGRI數(shù)據(jù)的反演結(jié)果符合實際的空氣污染發(fā)展趨勢及時空分布情況。

3.2 氣溶膠光學(xué)厚度產(chǎn)品交叉驗證

為了驗證研究方法在山東地區(qū)的適用性，針對1月22—24日反演的3個時次的結(jié)果，將其與JAXA官方發(fā)布的臨近時次AOD產(chǎn)品進(jìn)行比較分析。圖3a、b、c分別為22日12：00、23日12：00、24日12：00的Himawari-8/AHI氣溶膠分布圖。

圖3 與FY-4A/AGRI對應(yīng)時次的Himawari-8/AHI氣溶膠光學(xué)厚度(色階)產(chǎn)品圖(a. 22日， b. 23日， c. 24日)Fig.3 Himawari-8/AHI AOD product (color scale) at the same time with FY-4A/AGRI product on 22 (a), 23 (b), and 24 (c) January 2021

由圖3可以看出，Himawari-8/AHI的氣溶膠光學(xué)厚度產(chǎn)品與FY-4A/AGRI的反演結(jié)果在空間上具有很高的一致性，都主要分布在魯西南以及部分魯中地區(qū)；并且在AOD高低值的分布上也基本一致，高值主要分布在魯西南地區(qū)，也都存在著無法反演的空白區(qū)域，在圖上呈現(xiàn)出白色。由圖3b可以較為明顯地看出，23日相對于其他兩日AOD值都較大，主要的高值區(qū)體現(xiàn)在魯西南地區(qū)，在圖上表現(xiàn)為大片紅色區(qū)域，且高值達(dá)到0.6及以上，與本研究的反演結(jié)果一致。與前一天相比，23日空氣污染現(xiàn)象比較明顯，這在本算法的反演結(jié)果和Himawari-8/AHI中均有體現(xiàn)。總體來說，兩組數(shù)據(jù)具有較高的一致性，空間上AOD高低值區(qū)域分布相似且具有相同的變化趨勢。

選取FY-4A/AGRI的反演結(jié)果，在FY-4A/AGRI和Himawari-8/AHI同時存在有效值的情況下，提取對應(yīng)的Himawari-8/AHI產(chǎn)品數(shù)據(jù)，將二者的氣溶膠光學(xué)厚度數(shù)據(jù)進(jìn)行線性回歸分析，統(tǒng)計結(jié)果如圖4所示。

圖4 2021年1月22—24日FY-4A/AGRI反演值與Himawari-8/AHI產(chǎn)品值對比散點圖(a. 22日， b. 23日， c. 24日；紅線和黑線分別表示線性擬合回歸線以及1:1線，R代表相關(guān)系數(shù)，N代表有效像元數(shù)，RMSE(root mean squared error)代表均方根誤差，下同)Fig.4 Scatter plot of comparison between FY-4A/AGRI retrieval values and Himawari-8/AHI product values on 22 (a), 23 (b), and 24 (c) January 2021 (red and black lines for fitted regression line and 1:1 line, respectively, R for relative coefficient, N for effective pixels, RMSE for root mean squared error, the same hereafter)

圖4a為22日12：00的散點圖，有效值點共1 375個，線性擬合方程為y=0.37x+0.81，相關(guān)性較高，相關(guān)系數(shù)為0.92，誤差較小。圖4b有效值點數(shù)為32 246個，相關(guān)系數(shù)為0.92。圖4c為1月24日12：00的散點圖，有效值點共502個，線性擬合方程為y=1.1x+0.19，相關(guān)性較高，R為0.91。

由FY-4A/AGRI反演的AOD 結(jié)果與Himawari-8/AHI氣溶膠產(chǎn)品的對比散點圖(圖4)可以看出，自定義氣溶膠模型的反演結(jié)果與Himawari-8/AHI的氣溶膠產(chǎn)品之間具有明顯的線性相關(guān)關(guān)系，相關(guān)系數(shù)都在0.9以上，反演結(jié)果可以反映山東省氣溶膠的空間分布情況及變化趨勢。

為進(jìn)一步驗證本算法的可行性，同樣選取FY-4A/AGRI的反演結(jié)果和對應(yīng)時刻的MODIS氣溶膠光學(xué)厚度產(chǎn)品進(jìn)行對比驗證分析。圖5a、b、c分別為2021年1月22、23、24日FY-4A/AGRI數(shù)據(jù)12：00反演值與臨近時刻的MODIS產(chǎn)品值的對比散點圖。圖5a對比了兩幅圖中對應(yīng)的1 510個點，相關(guān)系數(shù)為0.93。圖5b中選取兩幅圖中的51 237個點，線性擬合方程為y=0.9x+0.11，相關(guān)系數(shù)為0.94。圖5c中有效值點共473個，線性擬合方程為y=0.71x+0.14，相關(guān)系數(shù)為0.94。相關(guān)系數(shù)也都在0.9以上，F(xiàn)Y-4A/AGRI數(shù)據(jù)反演值與MODIS產(chǎn)品具有顯著相關(guān)性。

圖5 2021年1月22—24日FY-4A/AGRI反演值與MODIS產(chǎn)品值對比散點圖(a. 22日， b. 23日， c. 24日)Fig.5 Scatter plot of comparison between FY-4A retrieval values and MODIS product values on 22 (a), 23 (b), and 24 (c) January 2021

圖6 2021年10月28—30日FY-4A反演AOD(色階)結(jié)果分布(a. 28日， b. 29日， c. 30日)Fig.6 Distribution of FY-4A retrieving AOD results (color scale) on 28 (a), 29 (b), and 30 (c) October 2021

針對山東省2021年1月的此次研究，F(xiàn)Y-4A/AGRI的氣溶膠反演結(jié)果與Himawari-8/AHI以及MODIS的氣溶膠產(chǎn)品均進(jìn)行對比驗證。對比結(jié)果表明，反演結(jié)果與兩者都具有較強(qiáng)的相關(guān)性，差異可能是因為基于FY-4A/AGRI進(jìn)行氣溶膠反演時，針對山東地區(qū)當(dāng)?shù)氐牡乩硖卣?、大氣條件等設(shè)置了更貼合實際情況的參數(shù)，自定義了山東地區(qū)的氣溶膠模型，而用來驗證的氣溶膠產(chǎn)品數(shù)據(jù)并沒有針對特定區(qū)域重新設(shè)定參數(shù)。

3.3 FY-4A衛(wèi)星2021年10月28—30日反演結(jié)果

通過查閱歷史資料可知，28日污染物開始逐漸累積，沿海地區(qū)以優(yōu)良天氣為主，魯西南地區(qū)以輕至中度污染為主，個別城市出現(xiàn)小時重度污染。29日山東省地面處于低壓場，擴(kuò)散條件較差，半島及魯南沿海地區(qū)以良至輕度污染為主，其他地區(qū)以輕至中度污染為主，其中魯西南、魯西北地區(qū)個別城市達(dá)到重度污染。30日地面處于弱氣壓場，氣流有輻合，不利于擴(kuò)散，污染物繼續(xù)累積，全省大部分地區(qū)以輕至中度污染為主，其中半島部分地區(qū)達(dá)到重度污染。

圖6分別為10月28日12：00、29日12：00、30日12：00的FY-4A/AGRI數(shù)據(jù)的氣溶膠光學(xué)厚度反演結(jié)果。10月28日AOD高值主要分布在魯西南地區(qū)，AOD較大值范圍為1.0～2.0，膠東半島區(qū)域AOD值較低。由圖6b可以看出，29日污染較為嚴(yán)重，整體的AOD值明顯增大，高值區(qū)域主要分布在魯西南和魯北地區(qū)，魯中區(qū)域的AOD值也呈現(xiàn)出增大的趨勢，在圖上顯示為大片的黃色高值區(qū)，空氣狀況進(jìn)一步惡化；與前一日相比，膠東半島地區(qū)的氣溶膠水平仍然較低，在圖上表現(xiàn)為藍(lán)色，空氣狀況比較樂觀。圖6c表明，30日除膠東半島地區(qū)，大部分區(qū)域都處于一個低值狀態(tài)，AOD值總體回歸較低水平，內(nèi)陸地區(qū)空氣污染狀況開始緩和。針對此次大氣污染的反演質(zhì)量較高，可以體現(xiàn)此次過程中AOD的分布情況以及變化趨勢，符合當(dāng)時的實際空氣質(zhì)量狀況。

3.4 Himawari-8/AHI氣溶膠光學(xué)厚度產(chǎn)品驗證

根據(jù)FY-4A/AGRI反演的結(jié)果顯示，在10月28—30日的空氣污染過程中，受影響區(qū)域最初在魯西南地區(qū)，后來影響范圍擴(kuò)大至魯西南及魯西北地區(qū)，最后高值區(qū)主要分布在膠東半島。為了驗證此次反演的可靠性，選取與AGRI反演數(shù)據(jù)對應(yīng)時刻的Himawari-8/AHI氣溶膠產(chǎn)品進(jìn)行對比分析，在空間分布上進(jìn)行直觀的比較。圖7分別為Himawari-8/AHI的10月28日12：00、29日12：00、30日12：00氣溶膠產(chǎn)品圖。

圖7 2021年10月28—30日與FY-4A/AGRI對應(yīng)時次的Himawari-8/AHI氣溶膠光學(xué)厚度(色階)產(chǎn)品圖(a. 28日， b. 29日， c. 30日)Fig.7 Himawari-8/AHI AOD product (color scale) at the same time with FY-4A/AGRI product on 28 (a), 29 (b), and 30 (c) October 2021

在圖7a中，28日的氣溶膠光學(xué)厚度高值區(qū)域同樣集中在魯西南地區(qū)，圖7b高值區(qū)域的分布與FY-4A/AGRI反演的結(jié)果一致，都集中在魯西北、魯西南地區(qū)。在圖7c中顯示30日的氣溶膠光學(xué)厚度分布情況也與FY-4A/AGRI反演結(jié)果相同，紅色的高值區(qū)域主要集中在膠東半島且部分黃色的高值區(qū)域在山東省的北部也有零星分布。

綜上所述，此次FY-4A/AGRI氣溶膠反演結(jié)果與Himawari-8/AHI的氣溶膠產(chǎn)品在空間分布上基本一致，并且具有相似的變化趨勢，為定量分析反演結(jié)果的準(zhǔn)確度，提取臨近時刻(分別為10月28日12：20、29日11：30、30日12：10)的MODIS產(chǎn)品數(shù)據(jù)，在保證FY-4A/AGRI和MODIS數(shù)據(jù)同時存在有效值的前提下，將兩者的氣溶膠光學(xué)厚度數(shù)據(jù)進(jìn)行線性回歸分析，統(tǒng)計結(jié)果如圖8所示。

在圖8a中，對比兩幅圖中對應(yīng)的19 159個點，線性擬合方程為y=0.85x+0.13，相關(guān)系數(shù)為0.9，具有較高的相關(guān)性。圖8b對比131 220個點，線性擬合方程為y=1.33x+0.35，相關(guān)系數(shù)為0.83。圖8c選取134 994個點，線性擬合方程為y=0.86x+0.46，相關(guān)系數(shù)為0.88。在研究過程中雖然盡最大可能選取時刻相近的數(shù)據(jù)進(jìn)行對比驗證分析，但二者之間仍然會有一定的時間差，而受濕度、溫度、風(fēng)速等氣象要素的影響，氣溶膠處于一種動態(tài)變化的狀態(tài)，這可能是導(dǎo)致結(jié)果存在差異的原因之一?？偟脕碚f，反演結(jié)果與MODIS氣溶膠產(chǎn)品的相關(guān)性較強(qiáng)，相關(guān)系數(shù)都在0.90左右。因此，可以認(rèn)為本研究方法具有一定的可靠性。

圖8 2021年10月28—30日FY-4A/AGRI反演值與MODIS產(chǎn)品值對比散點圖(a. 28日， b. 29日， c. 30日)Fig.8 Scatter plot of comparison between FY-4A/AGRI retrieval values and MODIS product values on 28 (a), 29 (b), and 30 (c) October 2021

由圖4、5、8可以看出，不同日期的反演結(jié)果分布存在差異，這可能是因為本研究的地表估算來自于固定的地表反射率數(shù)據(jù)庫，而葵花衛(wèi)星是多通道組合算法，二者估算出的結(jié)果本就存在差異，所以反演結(jié)果也會有所不同，但由于葵花衛(wèi)星沒有提供地表反射率結(jié)果，因而無法對這種分布差異進(jìn)行分析。針對AGRI的反演結(jié)果，選取對應(yīng)時刻對應(yīng)區(qū)域的MODIS氣溶膠產(chǎn)品數(shù)據(jù)，通過對比驗證兩組數(shù)據(jù)的AOD結(jié)果(表3)，發(fā)現(xiàn)相關(guān)系數(shù)均在0.9左右，具有較好的相關(guān)性。其中有4組數(shù)據(jù)進(jìn)行對比的有效像元數(shù)在1萬以上，分別是2021年1月23日，10月28、29、30日的數(shù)據(jù)，其余兩組數(shù)據(jù)的有效像元數(shù)量較少，可能是當(dāng)天的云較多，通過云檢測去除較多無法反演的像元，或是當(dāng)天的數(shù)據(jù)存在缺測等原因。對這兩組數(shù)據(jù)進(jìn)行統(tǒng)計分析發(fā)現(xiàn)，基于FY-4A/AGRI的自定義氣溶膠模型的反演結(jié)果大部分情況下低于MODIS的AOD值，根據(jù)已有的研究來看，可能是因為MODIS的AOD值存在高估的情況。

不論是由空間分布、變化趨勢的定性研究還是定量的線性回歸分析來看，F(xiàn)Y-4A/AGRI衛(wèi)星的AOD反演結(jié)果和氣溶膠產(chǎn)品都具有相似性和一致性，反演結(jié)果可以反映山東省氣溶膠光學(xué)厚度值的分布情況，并據(jù)此來對山東省的大氣環(huán)境污染進(jìn)行分析，監(jiān)測氣溶膠的時空分布。

表3 兩次天氣過程FY-4A/AGRI 的反演結(jié)果與MODIS氣溶膠產(chǎn)品對比結(jié)果

4 小結(jié)

以山東為研究區(qū)域，基于FY-4A/AGRI衛(wèi)星資料，充分考慮其通道設(shè)置，自定義山東地區(qū)的氣溶膠模型，研究其氣溶膠光學(xué)厚度反演算法，并將反演結(jié)果與臨近時刻的AOD產(chǎn)品進(jìn)行精度評價，主要結(jié)果如下：

(1)AGRI數(shù)據(jù)的反演結(jié)果與臨近時刻的Himawari-8/AHI產(chǎn)品在空間分布及變化趨勢上總體一致，且連續(xù)性、空間細(xì)節(jié)表達(dá)較好，可以用來監(jiān)測一天當(dāng)中AOD的動態(tài)變化情況。

(2)以MODIS氣溶膠產(chǎn)品為參考開展反演結(jié)果的定量分析研究，結(jié)果發(fā)現(xiàn)，AGRI的反演結(jié)果與臨近時刻的MODIS氣溶膠產(chǎn)品具有較好的相關(guān)性，相關(guān)系數(shù)均在0.90左右，說明本文的反演結(jié)果具有一定的合理性和可靠性，可以為大氣環(huán)境規(guī)劃與調(diào)控、污染防治提供決策依據(jù)。

(3)基于新一代國產(chǎn)靜止衛(wèi)星構(gòu)建氣溶膠反演算法，具有高時間分辨率的特點，能夠?qū)崿F(xiàn)高頻次實時動態(tài)的氣溶膠監(jiān)測，以此次研究成果為基礎(chǔ)，未來可以開展山東地區(qū)大氣氣溶膠的日變化特征以及大氣污染的時空變化監(jiān)測等方面的研究。