資源一號04星WFI數(shù)據(jù)氣溶膠反演與驗證
——以北京市及北京周邊地區(qū)為例

丁宇，汪小欽，王峰

(1. 數(shù)字中國研究院(福建)，福州 350108；2.衛(wèi)星空間信息技術(shù)綜合應(yīng)用國家地方聯(lián)合工程研究中心，福州 350108；3.福州大學(xué) 空間數(shù)據(jù)挖掘和信息共享教育部重點實驗室，福州 350108)

0 引言

氣溶膠是大氣層的重要組成成分，對地氣系統(tǒng)、輻射平衡、氣候變化有著顯著影響[1]，高濃度的氣溶膠也會加大人體患病的幾率[2]。傳統(tǒng)的氣溶膠觀測手段主要以地面實時監(jiān)測為主，難以滿足大范圍監(jiān)測和時空分布的要求，而遙感監(jiān)測手段彌補了這一不足，能夠?qū)Υ髿馕廴具M行大范圍的動態(tài)監(jiān)測[3]。

大氣氣溶膠遙感實質(zhì)上是研究如何實現(xiàn)地氣解耦，而實現(xiàn)地氣解耦的前提是得到準確的地表反射率信息。為此，研究者提出了暗像元[4]、偏振[5]、雙星協(xié)同[6]、深藍[7-8]等算法計算地表反射率，開展氣溶膠反演研究。其中，暗像元法和深藍算法在國內(nèi)外MODIS[9-10]、OLI[11-12]、HJ-1[13-14]、GF-1[15-16]、FY-3[17]、Himawar-8[18]等傳感器上得到廣泛應(yīng)用。暗像元法只對濃密植被等暗地表效果良好，深藍算法主要針對城市、沙漠等高亮地表。為了解決不能同時反演包含明暗像元的復(fù)雜地表上空氣溶膠光學(xué)厚度(aerosol optical depth，AOD)的問題，美國國家航空天局(National Aeronautics and Space Administration，NASA)最新發(fā)布的MODIS氣溶膠產(chǎn)品C6版本包括了融合2種算法的產(chǎn)品。文獻[19-20]將其與單一算法的MODIS氣溶膠產(chǎn)品對比，發(fā)現(xiàn)融合后的產(chǎn)品更具空間完整性。國內(nèi)學(xué)者結(jié)合2種算法在局部區(qū)域開展氣溶膠反演研究[21-23]，也得到了相似的結(jié)論。

資源一號04星(CBERS04)由我國與巴西聯(lián)合研制，并于2014年12月7日在太原衛(wèi)星發(fā)射中心發(fā)射成功。其搭載的寬視場成像儀WFI有B1(0.45～0.52 μm)、B2(0.52～0.59 μm)、B3(0.63～0.69 μm)和B4(0.77～0.89 μm)共4個波段，空間分辨率為73 m，幅寬為866 km，重訪周期高達3 d[24]。CBERS04WFI的時間分辨率和幅寬均優(yōu)于HJ-1 CCD、GF-1 WFV、OLI等傳感器，空間分辨率又優(yōu)于MODIS、FY-3 MERSI、Himawar-8 AHI等傳感器，但是目前利用WFI數(shù)據(jù)的應(yīng)用較少，還未見到應(yīng)用CBERS04WFI數(shù)據(jù)開展AOD反演的公開報道。

WFI在時空分辨率和幅寬上的優(yōu)勢，使其在氣溶膠估算上具有較大的應(yīng)用潛力。本文結(jié)合暗像元法與深藍算法，探討基于動態(tài)查找表的WFI數(shù)據(jù)AOD反演方法研究，為區(qū)域大氣變化監(jiān)測提供較高時空間分辨率的數(shù)據(jù)支持，同時擴展國產(chǎn)衛(wèi)星的應(yīng)用領(lǐng)域，拓展應(yīng)用廣度。

1 研究區(qū)和數(shù)據(jù)

1.1 研究區(qū)

本文以北京市及周邊地區(qū)(115°25′E～117°30′E，39°26′N～41°03′N)作為研究區(qū)。研究區(qū)位于華北平原西北部，研究區(qū)東南為平原，西北為山地，總體地勢為西北高、東南低(圖1)。研究區(qū)氣候?qū)倥瘻貛О霛駶櫚敫珊导撅L(fēng)氣候，年降水量時空分布不均勻，季節(jié)上夏多冬少，東北部和西南部山前迎風(fēng)坡地區(qū)為相對降水中心，西北部和北部深山區(qū)降水較少[25]。因此一年中北京6—8月份的大氣污染物的稀釋和擴散較快，大氣環(huán)境質(zhì)量相對最好；冬季受氣象條件和采暖的影響，春季受西北部沙塵暴的影響，大氣環(huán)境質(zhì)量相對較差[26]。

圖1 研究區(qū)海拔及AERONET站點分布

1.2 數(shù)據(jù)

研究所用的數(shù)據(jù)包括CBERS04WFI數(shù)據(jù)、MODIS地表反射率數(shù)據(jù)和氣溶膠產(chǎn)品，以及AERONET地基站點觀測數(shù)據(jù)。

從中國資源衛(wèi)星應(yīng)用中心網(wǎng)站下載的2016—2018年113景CBERS04星WFI數(shù)據(jù)，包括tif文件、觀測幾何txt文件，和提供定標參數(shù)等信息的xml文件。

2016—2018年的MODIS產(chǎn)品，包括MOD09A1數(shù)據(jù)和MODIS氣溶膠產(chǎn)品。MOD09A1數(shù)據(jù)是8 d合成的地表反射率產(chǎn)品，盡量降低了反射率短期變化以及云覆蓋的影響[27]。MODIS氣溶膠產(chǎn)品包括MOD04_3K(暗像元法，空間分辨率3 km)和MOD04_L2(深藍+暗像元，空間分辨率10 km)，為獲得較高時空分辨率的氣溶膠產(chǎn)品，將2種產(chǎn)品進行融合，融合原則是以MOD04_3K數(shù)據(jù)為主，用MOD04_L2數(shù)據(jù)補充缺失區(qū)域。

采用AERONET站點Level 2.0級數(shù)據(jù)驗證反演結(jié)果。因AERONET數(shù)據(jù)沒有與反演結(jié)果對應(yīng)的550 nm處AOD數(shù)據(jù)，利用二次多項式插值得到550 nm處AOD值作為驗證數(shù)據(jù)[28]，如式(1)所示。

lnτλ=τλ+a1lnλ+a2(lnλ)2

(1)

式中：τλ表示λnm處的AOD值；a1、a2為相應(yīng)項未知系數(shù)。

2 研究方法

2.1 基本原理

假設(shè)地面為均勻的朗伯面且大氣條件均一，衛(wèi)星傳感器的表觀反射率為ρTOA，如式(2)所示。

(2)

式中：ρO為大氣的路徑輻射項等效反射率；ρs表示地表二向反射率；φ為相對方位角；μs=cosθs，μv=cosθv，μs與μv為太陽天頂角和觀測天頂角；S為大氣下界的半球反射率；T(μs)、T(μv)分別為大氣分子上下行散射透過率[29]。

通過式(2)可以明顯看出，要想獲得準確的大氣路徑輻射項等效反射率，需要計算得到準確的地表反射率。因此，本研究結(jié)合深藍算法和暗像元法，對亮暗像元進行區(qū)分，分別計算其地表反射率值。其中，暗像元法是根據(jù)濃密植被等暗地物在紅藍波段與短波紅外波段的地表反射率間的線性關(guān)系實現(xiàn)氣溶膠反演。但是由于眾多傳感器不具有短波紅外波段，無法使用該算法進行AOD反演。因此許多研究基于紅藍通道固定地表反射率比值法進行氣溶膠反演，式(3)為NASA的MODIS 數(shù)據(jù)V5.2氣溶膠算法的紅藍通道地表反射率關(guān)系式[30]。

y=0.49x+0.05

(3)

式中：y為紅光波段反射率；x為藍光波段反射率。因WFI傳感器缺失短波紅外波段信息，選擇利用紅藍波段間線性關(guān)系開展WFI傳感器氣溶膠反演研究。為得到WFI傳感器紅藍波段間線性關(guān)系，將WFI和MODIS紅藍波段的光譜響應(yīng)函數(shù)與ENVI光譜庫中220種植被光譜(重采樣為1 nm)進行卷積，計算WFI和MODIS紅藍通道植被地表反射率并進行回歸分析(圖2)。結(jié)合圖2中擬合方程式和式(3)，計算得到WFI傳感器紅藍波段地表反射率關(guān)系模型，如式(4)所示。

圖2 MODIS與WFI紅藍波段地表反射率散點圖

ρ3=0.500 4ρ1+0.046 8

(4)

式中：ρ1和ρ3分別為WFI藍光與紅光地表反射率。

針對高亮的城市、沙漠等區(qū)域，本研究選擇采用的深藍算法是依據(jù)藍光波段氣溶膠對衛(wèi)星觀測信號有較強貢獻而地表反射較弱的特點，借助地表反射率庫先驗知識，以實現(xiàn)地氣解耦。本文選擇使用MOD09A1數(shù)據(jù)來構(gòu)建地表反射率庫，但不同傳感器對相同地物成像得到的地表反射率會有差異。研究表明，0.01的地表反射誤差帶來約0.1的氣溶膠光學(xué)厚度誤差[31]。因此，需要通過波段修正的方法將MOD09A1數(shù)據(jù)轉(zhuǎn)換為WFI地表反射率數(shù)據(jù)。選用ENVI光譜庫中各類地物光譜共426條，分別與MODIS和WFI藍光光譜響應(yīng)函數(shù)卷積，計算光譜等效值并進行回歸分析(圖3)。利用二者擬合方程式將MOD09A1數(shù)據(jù)轉(zhuǎn)換為WFI地表反射率數(shù)據(jù)。如圖3所示，光譜修正公式與直線y=x極其接近，因此轉(zhuǎn)換得到的WFI地表反射率與經(jīng)過質(zhì)量驗證的MOD09A1地表反射率產(chǎn)品只有極小的差異，但經(jīng)過光譜修正得到的WFI地表反射率更符合WFI傳感器光譜特性。

圖3 MODIS與WFI藍波段地表反射率擬合圖

2.2 WFI數(shù)據(jù)氣溶膠反演流程

遙感反演氣溶膠光學(xué)厚度的主要過程是將相關(guān)先驗知識，如觀測幾何、氣溶膠模型、觀測時間等，迭代循環(huán)輸入6SV輻射傳輸模型，計算得到3個大氣參數(shù)ρo、S和T。結(jié)合暗像元法或深藍算法計算得出的地表反射率參與輻射傳輸計算，得到最終AOD值，技術(shù)路線如圖4所示。

圖4 CBERS04星WFI數(shù)據(jù)氣溶膠反演流程圖

1)輻射定標。利用原始數(shù)據(jù)中的xml文件自帶的輻射定標系數(shù)、太陽高度角和成像時間等參數(shù)，結(jié)合大氣層外太陽輻照度[32]，將影像DN值轉(zhuǎn)換為表觀反射率，如式(5)所示。

(5)

式中：ρ為表觀反射率；Gain為增益值；Bias為偏置值；D為日地距離；ESUNb為大氣層頂太陽輻照度；θ為太陽天頂角[33]。

2)云水檢測。影像上的云、水像元會對反演結(jié)果造成誤差，需要在實驗前對云、水進行檢測。將歸一化水體指數(shù)NDWI大于0.3的像元判定為水體像元[34]，將紅光波段大于0.25且歸一化植被指數(shù)NDVI[35]<0的像元判斷為云像元[23]。

3)亮暗像元判斷。通過增強型植被指數(shù)EVI[36]區(qū)分亮暗像元。將EVI大于0.4的像元判定為暗像元[37]，暗像元區(qū)域采用暗像元法反演AOD；EVI小于等于0.4的亮像元區(qū)域則采用深藍算法(式(6))。

(6)

式中：ρ1為藍光波段值；ρ3為紅光波段值；ρ4為近紅外波段值。

4)動態(tài)查找表構(gòu)建。傳統(tǒng)查找表只針對單個傳感器，其構(gòu)建方法是事先設(shè)置好參數(shù)組合，通過大氣輻射傳輸模型計算得到，可移植性不強。動態(tài)查找表是直接在反演過程中，根據(jù)每一景影像的幾何參數(shù)生成只針對該影像的查找表。查找表構(gòu)建成功后，以臨時文件的形式保存在計算機內(nèi)存中，查找效率遠超以文本形式保存的傳統(tǒng)查找表。動態(tài)查找表可以輸入多個傳感器的光譜響應(yīng)函數(shù)，同時適用多個傳感器，具有較強的可移植性。構(gòu)建方法：①將影像太陽、衛(wèi)星天頂角數(shù)據(jù)四舍五入為整數(shù)值(WFI傳感器的觀測幾何數(shù)據(jù)中這2個角度覆蓋范圍最大)，按照角度遞增順序劃分區(qū)域；②對2個角度文件劃分區(qū)域取并集，得到最終劃分區(qū)域；③逐個區(qū)域計算觀測幾何參數(shù)，將預(yù)設(shè)AOD值(550 nm處AOD步長設(shè)置范圍為0～1.95之間，步長為0.05)，以及氣溶膠模型輸入6SV輻射傳輸模型計算得到大氣參數(shù)數(shù)據(jù)。本文選擇能夠較準確描述研究區(qū)氣溶膠微粒狀況的大陸型氣溶膠模型，作為6SV模型輸入條件，該模型提供了所需的散射相函數(shù)、單次散射反照率和不對稱因子等關(guān)鍵參數(shù)[38]。

3 結(jié)果與分析

按照上述算法與數(shù)據(jù)處理流程，反演得到研究區(qū)空間分辨率為500 m×500 m的AOD分布。

3.1 氣溶膠反演方法對比與分析

為分析本文方法應(yīng)用于反演復(fù)雜地表上空AOD的性能，分別與暗像元法和深藍算法結(jié)果進行比較(圖5)。從圖中可以發(fā)現(xiàn)：暗像元法較好地反演了濃密植被覆蓋的山區(qū)，但無法反演亮像元的城市等區(qū)域；深藍算法能夠有效反演亮地表上空AOD，但對于其他區(qū)域的反演結(jié)果偏高，其反演結(jié)果均值為0.223 1，而本文方法反演結(jié)果均值為0.177 7，說明深藍算法反演結(jié)果整體偏高，這也與張勝敏等[39]的研究結(jié)論一致。究其原因是因為深藍算法使用的MOD09A1地表反射率產(chǎn)品是8 d數(shù)據(jù)合成構(gòu)建而成的，這一構(gòu)建方法會低估地表反射率，進而在反演中高估大氣反射率，最終造成反演結(jié)果偏高[27]。本文方法有效集成了暗像元法和深藍算法的優(yōu)勢，反演結(jié)果較單一算法更適合用于反演包含明暗像元的復(fù)雜地表上空AOD。

圖5 2018年9月21日的反演結(jié)果

3.2 與MODIS產(chǎn)品對比

為較為全面地驗證反演精度，本研究使用融合后的MODIS氣溶膠產(chǎn)品與反演結(jié)果進行空間分布的對比。如圖6所示，反演得到的 AOD與MODIS氣溶膠產(chǎn)品高、低值分布情況基本一致，但覆蓋范圍更全；MOD04_3K(暗像元法)和MOD04_L2(深藍+暗像元)氣溶膠產(chǎn)品，都會在生產(chǎn)過程中去除不合適的像元(如云、沙漠、雪/冰和內(nèi)陸水)，同時也剔除一部分最暗和最亮像元[40]，因此，融合后的MODIS產(chǎn)品還是會存在較多無效值，在覆蓋范圍上不如WFI AOD；AOD值最高的是研究區(qū)東南部的通州、大興、廊坊等地，中部區(qū)域次之，西北部房山、門頭溝、昌平、延慶、懷柔、密云等區(qū)域AOD值最低，即研究區(qū)AOD分布總體上呈現(xiàn)自東南城市區(qū)域向西北山區(qū)遞減的趨勢。與圖1對比發(fā)現(xiàn)，AOD分布情況與海拔存在明顯的對應(yīng)關(guān)系，人類活動多的低海拔平坦區(qū)域AOD較高，山區(qū)AOD較低。

圖6 WFI反演AOD與MODIS AOD空間分布對比圖

3.3 基于AERONET站點監(jiān)測數(shù)據(jù)的精度驗證

本文采用AERONET站點Level 2.0級數(shù)據(jù)驗證反演結(jié)果。為客觀評價研究模型反演結(jié)果的精度，本文采用了相關(guān)系數(shù)(r)、均方根誤差(RMSE)、平均絕對誤差(MAE)和預(yù)期誤差區(qū)間(expected error，EE)這4個統(tǒng)計指標對其進行精度評定。其中，預(yù)期誤差區(qū)間EE的計算公式如式(7)所示[41]。

EE=±(0.05+0.2·τα)

(7)

式中：τα為AERONET站點實測AOD值。

如圖7所示，本研究共得220個有效數(shù)據(jù)對，反演結(jié)果與AERONET站點實測值相關(guān)性較高(r>0.94)，RMSE為0.170 3，MAE為0.118 3，67%的反演結(jié)果在誤差區(qū)間內(nèi)，29%的反演結(jié)果高估，4%的反演結(jié)果低估。

圖7 反演結(jié)果與AERONET AOD散點對比圖

高估現(xiàn)象主要發(fā)生在低值A(chǔ)OD處，即AERONET站點AOD小于0.4時。這一現(xiàn)象主要有以下2點原因造成：①地表反射率相差0.01，AOD會出現(xiàn)約0.1的偏差[31]。因此在實際AOD值較低時，地表反射率微弱的偏差都可能導(dǎo)致反演結(jié)果較大的偏差。②高估的原因可能與MODIS地表反射率數(shù)據(jù)構(gòu)建方式有關(guān)。MOD09A1地表反射率數(shù)據(jù)的構(gòu)建盡量考慮高觀測覆蓋、低視角、無云及云的陰影以及氣溶膠濃度的影響，選擇8 d中最可能、最穩(wěn)定的數(shù)據(jù)構(gòu)建而成[27]。而在大氣質(zhì)量較好時，通過這種方式構(gòu)建出的地表反射率數(shù)據(jù)，可能會因為考慮太多因素導(dǎo)致反射率略微偏低，從而導(dǎo)致氣溶膠反演結(jié)果偏高。

綜上所述，利用本文方法反演氣溶膠光學(xué)厚度具有較高穩(wěn)定性，反演結(jié)果的空間分布和精度都符合要求。

4 結(jié)束語

本文基于CBERS04星WFI數(shù)據(jù)特點，結(jié)合暗像元法和深藍算法，實現(xiàn)了基于動態(tài)查找表的AOD估算。本文利用WFI數(shù)據(jù)結(jié)合暗像元法和深藍算法較好地反演了包含亮暗像元的復(fù)雜地表上空AOD，避免了單一算法的不足。本文反演結(jié)果與MODIS氣溶膠產(chǎn)品和AERONET站點AOD進行了對比驗證，證明本文研究方法能較好地反映氣溶膠的空間分布。反演結(jié)果空間分辨率和空間覆蓋都優(yōu)于MODIS氣溶膠產(chǎn)品，且與地面站點實測值相關(guān)性顯著(r>0.94)，但在空氣質(zhì)量較好時反演結(jié)果容易高估，這與地表反射率數(shù)據(jù)的精度及其構(gòu)建方式有關(guān)。