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

陳 瑞 游浩妍 萬 翔

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

0 引言

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

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

1 資料與方法

1.1 資料

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

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

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

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

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

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

1.2 反演原理

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

(1)

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

(2)

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

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

1.3 反演方法

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

1.3.1云掩膜

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

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

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

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

式中

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

(9)

式中,

(10)

1.3.3查找表建立

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

表1 查找表示例

1.3.4查找表插值

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

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

2 結(jié)果與分析

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

2.1 與MODIS結(jié)果對比

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

從回歸分析結(jié)果可見: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é)果分布對比

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

2.2 與AERONET站點實測結(jié)果對比

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

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

3 結(jié)束語

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