基于Himawari-8衛星的氣溶膠光學厚度反演

陳 瑞 游浩妍 萬 翔

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

0 引言

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

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

1 資料與方法

1.1 資料

1.1.1Himawari-8數據

H8數據來源于日本氣象廳,于2015年7月正式投入使用,用戶通過網站注冊后可以免費獲取數據(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,特定地區編碼掃描可達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數據

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

1.1.3AERONET數據

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

1.2 反演原理

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

(1)

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

(2)

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

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

1.3 反演方法

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

1.3.1云掩膜

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

式中,σ0.51為標準差;ρ0.51為H8數據波長0.51 μm波段;ρ0.47為H8數據波長0.47 μm波段;μ0.51為窗口內H8波長0.51 μm波段的平均值。

1.3.2地表波段關系

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

式中

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

(9)

式中,

(10)

1.3.3查找表建立

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

表1 查找表示例

1.3.4查找表插值

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

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

2 結果與分析

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

2.1 與MODIS結果對比

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

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

(a)UTC 03 ：50 H8反演

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

(c)UTC 05 ：30 H8反演

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

圖1 2020年5月1日H8反演AOD結果與MODIS AOD結果分布對比

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

2.2 與AERONET站點實測結果對比

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

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

3 結束語

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