基于Himawari-8衛(wèi)星的氣溶膠光學(xué)厚度反演

陳 瑞 游浩妍 萬(wàn) 翔

(1. 自然資源部第一地理信息制圖院, 陜西 西安 710054;2. 自然資源部陜西基礎(chǔ)地理信息中心, 陜西 西安 710054)

0 引言

氣溶膠是大氣系統(tǒng)重要構(gòu)成,對(duì)氣候變化、人類健康有重要影響,在全球和區(qū)域氣候變化中扮演十分重要的角色[1-2]。氣溶膠光學(xué)厚度是描述光的衰減作用的無(wú)量綱參數(shù),能表征大氣渾濁度,是大氣氣溶膠基本光學(xué)特性之一[3]。衛(wèi)星遙感可以開(kāi)展氣溶膠光學(xué)特性的宏觀監(jiān)測(cè),為獲取氣溶膠時(shí)空分布特征、變化趨勢(shì)及溯源信息提供有效監(jiān)測(cè)手段。

高時(shí)間分辨率是地球靜止衛(wèi)星的優(yōu)勢(shì),能實(shí)現(xiàn)分鐘級(jí)的對(duì)地觀測(cè),對(duì)氣溶膠的日變化情況研究具有巨大的應(yīng)用潛力[4-6]。葵花8號(hào) (Himawari-8,H8)是地球靜止衛(wèi)星的代表,于2015年7月開(kāi)始向全球?qū)崟r(shí)廣播數(shù)據(jù)[7]。目前基于H8數(shù)據(jù)開(kāi)展的氣溶膠監(jiān)測(cè)已有相關(guān)研究:文獻(xiàn)[8]利用暗像元法開(kāi)展了H8氣溶膠反演,并對(duì)結(jié)果進(jìn)行驗(yàn)證,精度可行;文獻(xiàn)[9]利用新研究算法對(duì)H8衛(wèi)星氣溶膠光學(xué)厚度進(jìn)行反演,并通過(guò)全球氣溶膠地基觀測(cè)網(wǎng)(aerosol robotic network，AERONET)站點(diǎn)進(jìn)行誤差分析,說(shuō)明通過(guò)H8反演的結(jié)果在霧霾傳輸、成因分析等方面有著廣泛應(yīng)用價(jià)值;文獻(xiàn)[10]構(gòu)建了基于最優(yōu)估計(jì)技術(shù)的H8衛(wèi)星氣溶膠反演算法,并驗(yàn)證了可行性。本文參考中分辨率成像光譜儀(moderate-resolution imaging spectroradiometer，MODIS)暗目標(biāo)法,構(gòu)建H8的氣溶膠光學(xué)厚度(aerosol optical depth，AOD)反演算法,并將結(jié)果與AERONET、MODIS AOD進(jìn)行對(duì)比,對(duì)反演結(jié)果進(jìn)行驗(yàn)證分析。

1 資料與方法

1.1 資料

1.1.1Himawari-8數(shù)據(jù)

H8數(shù)據(jù)來(lái)源于日本氣象廳,于2015年7月正式投入使用,用戶通過(guò)網(wǎng)站注冊(cè)后可以免費(fèi)獲取數(shù)據(jù)(https://www.eorc.jaxa.jp/ptree/)。搭載的高級(jí)成像儀(Advanced Himawari Imager,AHI)傳感器,包括16個(gè)通道,覆蓋可見(jiàn)光、近紅外、紅外波段,星下點(diǎn)空間分辨率第3波段為0.5 km,1、2、4波段為1 km,5～16波段為2 km[11]。先進(jìn)葵花成像儀(advanced himawari imager，AHI)全圓盤圖完成周期是10 min,特定地區(qū)編碼掃描可達(dá)2.5 min。本文利用H8的1～6波段進(jìn)行AOD反演,其中云檢測(cè)利用第1(0.46 μm)、2(0.51 μm)、4(0.86 μm)、5(1.6 μm)波段進(jìn)行計(jì)算;AOD利用1(0.46 μm)、3(0.64 μm)、6(2.3 μm)波段進(jìn)行反演。

1.1.2MOD04數(shù)據(jù)

MOD04數(shù)據(jù)是美國(guó)國(guó)家航空航天局(National Aeronautics and Space Administration，NASA)在國(guó)際上公開(kāi)發(fā)布的包含各類氣溶膠特性參數(shù)的AOD產(chǎn)品,數(shù)據(jù)在官網(wǎng)(https://ladsweb.modaps.eosdis.nasa.gov/search/)免費(fèi)獲得。產(chǎn)品經(jīng)過(guò)大量地面實(shí)測(cè)站點(diǎn)驗(yàn)證,具有較高的精度,空間分辨率3 km,一天過(guò)境2次,過(guò)境時(shí)間為地方時(shí)10 ：30和13 ：30左右,難以滿足小區(qū)域氣溶膠的日變化監(jiān)測(cè)分析,常作為對(duì)比產(chǎn)品來(lái)檢驗(yàn)反演結(jié)果的精度。本文利用MOD04產(chǎn)品中暗像元法反演的氣溶膠產(chǎn)品,通過(guò)時(shí)空匹配作為H8 AOD反演結(jié)果的對(duì)比驗(yàn)證數(shù)據(jù),空間分辨率為3 km。

1.1.3AERONET數(shù)據(jù)

AERONET可以采集更為精確的氣溶膠特性數(shù)據(jù),是NASA聯(lián)合其他科研機(jī)構(gòu)等設(shè)立的自動(dòng)觀測(cè)網(wǎng),提供了全球范圍內(nèi)多個(gè)站點(diǎn)的氣溶膠光學(xué)厚度數(shù)據(jù)[12],利用CE-318型太陽(yáng)輻射計(jì)以及統(tǒng)一的算法獲取AOD信息,數(shù)據(jù)可從官網(wǎng)免費(fèi)獲得(https://aeronet.gsfc.nasa.gov/)。本文利用北京區(qū)域AERONET長(zhǎng)期觀測(cè)站點(diǎn):Beijing、Beijing_CAMS、XiangHe,采用經(jīng)人工篩選質(zhì)量級(jí)別高的Level2.0產(chǎn)品作為驗(yàn)證數(shù)據(jù),利用440～675nm波長(zhǎng)處的AOD,插值得到550nm波段的AOD。

1.2 反演原理

衛(wèi)星獲取的表觀反射率變化能夠體現(xiàn)氣溶膠消光作用,從而進(jìn)行AOD反演。根據(jù)大氣輻射傳輸方程,在地表郎伯體和大氣水平均一的前提下,表觀反射率ρTOA可以表示為式(1)。

(1)

式中,μS,μV,φ分別為太陽(yáng)天頂角余弦、觀測(cè)天頂角余弦、相對(duì)方位角;ρ0是大氣中分子和氣溶膠散射的反射率;ρS是地表反射率；T(μS)是入射方向的大氣透過(guò)率；T(μV)是觀測(cè)方向的大氣透過(guò)率;S為大氣下界面的半球反射率。對(duì)于單次散射,ρ0與氣溶膠光學(xué)厚度τ、氣溶膠散射相函數(shù)Pa、單次散射反照率ω0以及分子散射造成的路徑輻射ρm有關(guān)[13],公式表示為式(2)。

(2)

由式(2)可見(jiàn),當(dāng)ρS很小時(shí),ρTOA主要取決于大氣貢獻(xiàn)項(xiàng);當(dāng)ρS很大時(shí),地表成為主要貢獻(xiàn)項(xiàng)。衛(wèi)星數(shù)據(jù)可以提供角度信息及表觀反射率信息,氣溶膠散射相函數(shù)Pa、單次散射反照率ω0以及路徑輻射ρm、大氣透過(guò)率T(μS)T(μV)、半球反射率S等參數(shù)我們通過(guò)6S模型計(jì)算獲取。因此,求得地表反射率后,衛(wèi)星觀測(cè)的表觀反射率除掉地表反射率部分即得到AOD值。

MODIS暗目標(biāo)法是在清潔大氣情況下,暗目標(biāo)地物(植被、水體、黑土等)在近紅外通道地表反射率很小,且近紅外波段與紅藍(lán)波段地表反射率具有較高的相關(guān)性,因此利用近紅外波段估算紅藍(lán)波段的地表反射率,通過(guò)6S模型建立查找表,實(shí)現(xiàn)AOD反演[14]。

1.3 反演方法

利用交互式編程語(yǔ)言實(shí)現(xiàn)AOD的反演,工藝如下:將紅外波段中數(shù)值小于0.25的像元選擇為暗像元,并得到紅藍(lán)波段相同位置的暗像元;根據(jù)MODIS C6版本中暗目標(biāo)法的地表反射率關(guān)系計(jì)算紅藍(lán)通道地表反射率;讀取每個(gè)像元角度信息,并尋找查找表中對(duì)應(yīng)的角度值上下限,對(duì)上下限角度進(jìn)行插值,得到大氣透過(guò)率T(μS)T(μV)、半球反射率S、路徑輻射ρm等大氣參數(shù);通過(guò)3個(gè)大氣參數(shù)、紅、藍(lán)波段地表反射率結(jié)合式(1),得出紅藍(lán)波段的理論表觀反射率。最接近實(shí)際表觀反射率的查找表AOD為反演結(jié)果,最后對(duì)紅藍(lán)波段的反演結(jié)果進(jìn)行平均。

1.3.1云掩膜

統(tǒng)計(jì)H8數(shù)據(jù)云覆蓋區(qū)域的表觀反射率情況,設(shè)置閾值進(jìn)行去云處理,當(dāng)3×3窗口的表觀反射率ρ0.47>0.25,或者表觀反射率ρ0.51標(biāo)準(zhǔn)差>0.006,則窗口內(nèi)的9個(gè)像元均認(rèn)為是云。

式中,σ0.51為標(biāo)準(zhǔn)差;ρ0.51為H8數(shù)據(jù)波長(zhǎng)0.51 μm波段;ρ0.47為H8數(shù)據(jù)波長(zhǎng)0.47 μm波段;μ0.51為窗口內(nèi)H8波長(zhǎng)0.51 μm波段的平均值。

1.3.2地表波段關(guān)系

MODIS C6版本中暗目標(biāo)法的地表反射率關(guān)系為式(5)和式(6)。

式中

H8波段設(shè)置同MODIS相似,利用MODIS波段間地表反射率關(guān)系進(jìn)行H8衛(wèi)星AOD反演。由于H8 AHI傳感器沒(méi)有1.24 μm通道,統(tǒng)計(jì)0.86 μm波段與1.24 μm波段關(guān)系,以替代1.24 μm計(jì)算的NDVI[15]。最終適用于H8的波段地表反射率關(guān)系見(jiàn)式(8)和式(9)。

(9)

式中,

(10)

1.3.3查找表建立

利用6S模型建立北京區(qū)域5月1日查找表,輸入?yún)?shù)包括:平均海拔0.47 km;太陽(yáng)天頂角為24°～68°,間隔1°;觀測(cè)天頂角為51°～55°,間隔1°;相對(duì)方位角為0°～180°,間隔12°;550 nmAOD為0～2.5,間隔0.1;選擇大陸型氣溶膠;根據(jù)上述參數(shù)以及6S模型模擬的大氣參數(shù)建立査找表,如表1所示。

表1 查找表示例

1.3.4查找表插值

查找表包含的角度值不完全等于像元讀取的角度,故對(duì)在查找表中查找到對(duì)應(yīng)角度的上下限進(jìn)行插值得到前面3個(gè)參數(shù)的插值結(jié)果,以此作為參數(shù)計(jì)算表觀反射率。分析查找表發(fā)現(xiàn):相同AOD的半球反照率相同,太陽(yáng)天頂角(Solar zenith angle, SOZ)和衛(wèi)星天頂角(satellite zenith angle，SAZ)同大氣總通過(guò)率成反比,與大氣中氣溶膠反射率成正比;相對(duì)方位角與大氣中氣溶膠反射率成反比,與大氣總透過(guò)率無(wú)關(guān)。插值方法見(jiàn)式(13)～式(16)。

式中,C為插值系數(shù);Z為太陽(yáng)天頂角;ZA為衛(wèi)星天頂角,φ為此像元讀取的角度數(shù)據(jù);S為半球反照率;T為大氣總透過(guò)率;P為大氣中氣溶膠反射率,通過(guò)查找表1得出,Z1、ZA1為查找出的角度數(shù)據(jù)上限,S1、T1、P1是上限對(duì)應(yīng)的半球反照率、大氣總透過(guò)率、大氣中氣溶膠反射率,Z2、ZA2、S2、T2、P2為查找出的角度數(shù)據(jù)下限。

2 結(jié)果與分析

本文以北京為研究區(qū)域,選擇2020年5月1日UTC時(shí)間00 ：00—08 ：00逐小時(shí)數(shù)據(jù)進(jìn)行反演,當(dāng)日植被情況較好、云量較少、有地面驗(yàn)證點(diǎn),可保證足夠的反演結(jié)果。當(dāng)日空氣質(zhì)量指數(shù)(air quality index，AQI)為204,屬于嚴(yán)重污染天氣。

2.1 與MODIS結(jié)果對(duì)比

H8衛(wèi)星反演AOD結(jié)果為2 km分辨率,重采樣為3 km分辨率后與MODIS暗像元的AOD產(chǎn)品結(jié)果進(jìn)行對(duì)比分析,從空間分布圖可見(jiàn):H8和MODIS的AOD結(jié)果空間分布一致,但存在細(xì)微差別：①來(lái)自云掩膜方法不同;②分辨率不同重采樣后進(jìn)行對(duì)比。

從回歸分析結(jié)果可見(jiàn):H8反演AOD結(jié)果與MODIS AOD結(jié)果相關(guān)系數(shù)為0.854,相關(guān)性較好,表明算法具有可行性,如圖1、圖2所示。

(a)UTC 03 ：50 H8反演

AOD結(jié)果 (b)UTC 03 ：45 MODIS AOD

(c)UTC 05 ：30 H8反演

AOD結(jié)果 (d)UTC 05 ：25 MODIS AOD

圖1 2020年5月1日H8反演AOD結(jié)果與MODIS AOD結(jié)果分布對(duì)比

圖2 2020年5月1日UTC 03 ：45及05 ：25 H8 AOD反演結(jié)果與MODIS AOD結(jié)果對(duì)比

2.2 與AERONET站點(diǎn)實(shí)測(cè)結(jié)果對(duì)比

圖3 2020年5月1日H8 AOD結(jié)果與AERONET站點(diǎn)對(duì)比

H8衛(wèi)星反演結(jié)果為550 nm處的AOD,而AERONET沒(méi)有直接觀測(cè)的數(shù)據(jù),根據(jù)其440 nm、500 nm和670 nm處的觀測(cè)結(jié)果,利用二次多項(xiàng)式擬合方法,計(jì)算得到550 nm處AOD直接觀測(cè)結(jié)果。結(jié)果對(duì)比時(shí),時(shí)間上將H8觀測(cè)時(shí)間前后5 min的地測(cè)結(jié)果平均值作為該時(shí)刻觀測(cè)值;空間上,以站點(diǎn)為中心,半徑10 km區(qū)域內(nèi)H8 AOD像元均值為該站點(diǎn)處AOD反演結(jié)果。驗(yàn)證結(jié)果如圖3所示,可見(jiàn)AOD反演結(jié)果和AERONET站點(diǎn)結(jié)果相關(guān)系數(shù)為0.809,相關(guān)性較高,與MODIS的對(duì)比分析情況相似,具有低值區(qū)AOD略高、高值區(qū)略低的趨勢(shì)。

3 結(jié)束語(yǔ)

本文構(gòu)建了H8的AOD反演算法,并以北京地區(qū)2020年5月1日為對(duì)象反演了AOD結(jié)果,通過(guò)與AERONET 3個(gè)站點(diǎn)Beijing-RADI、Beijing-CAMS、Xianghe的觀測(cè)數(shù)據(jù)(MODIS AOD產(chǎn)品)進(jìn)行對(duì)比,結(jié)果具有較好的相關(guān)性和驗(yàn)證精度,可見(jiàn)利用暗目標(biāo)法對(duì)H8衛(wèi)星反演AOD具有較大的應(yīng)用潛力。相比于MODIS AOD產(chǎn)品,該算法時(shí)間分辨率高,能得到10 min級(jí)、空間分辨率2 km的AOD日變化結(jié)果,對(duì)于動(dòng)態(tài)監(jiān)測(cè)氣溶膠情況具有較好的潛力。存在的問(wèn)題及下一步工作包括:①目前查找表構(gòu)建時(shí)采用大陸型氣溶膠,下一步通過(guò)研究區(qū)域氣溶膠實(shí)際構(gòu)成情況,采用自定義氣溶膠模式,提升AOD反演結(jié)果精度;②因地面觀測(cè)站點(diǎn)有限,將MODIS地表反射率關(guān)系轉(zhuǎn)化為H8衛(wèi)星,后續(xù)可進(jìn)一步優(yōu)化波段地表反射率關(guān)系,實(shí)現(xiàn)地表反射率的精確估計(jì)。