丁宇，汪小欽，王峰

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

0 引言

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

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

資源一號(hào)04星(CBERS04)由我國(guó)與巴西聯(lián)合研制，并于2014年12月7日在太原衛(wèi)星發(fā)射中心發(fā)射成功。其搭載的寬視場(chǎng)成像儀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個(gè)波段，空間分辨率為73 m，幅寬為866 km，重訪周期高達(dá)3 d[24]。CBERS04WFI的時(shí)間分辨率和幅寬均優(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反演的公開報(bào)道。

WFI在時(shí)空分辨率和幅寬上的優(yōu)勢(shì)，使其在氣溶膠估算上具有較大的應(yīng)用潛力。本文結(jié)合暗像元法與深藍(lán)算法，探討基于動(dòng)態(tài)查找表的WFI數(shù)據(jù)AOD反演方法研究，為區(qū)域大氣變化監(jiān)測(cè)提供較高時(shí)空間分辨率的數(shù)據(jù)支持，同時(shí)擴(kuò)展國(guó)產(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ū)東南為平原，西北為山地，總體地勢(shì)為西北高、東南低(圖1)。研究區(qū)氣候?qū)倥瘻貛О霛駶?rùn)半干旱季風(fēng)氣候，年降水量時(shí)空分布不均勻，季節(jié)上夏多冬少，東北部和西南部山前迎風(fēng)坡地區(qū)為相對(duì)降水中心，西北部和北部深山區(qū)降水較少[25]。因此一年中北京6—8月份的大氣污染物的稀釋和擴(kuò)散較快，大氣環(huán)境質(zhì)量相對(duì)最好；冬季受氣象條件和采暖的影響，春季受西北部沙塵暴的影響，大氣環(huán)境質(zhì)量相對(duì)較差[26]。

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

1.2 數(shù)據(jù)

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

從中國(guó)資源衛(wèi)星應(yīng)用中心網(wǎng)站下載的2016—2018年113景CBERS04星WFI數(shù)據(jù)，包括tif文件、觀測(cè)幾何txt文件，和提供定標(biāo)參數(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(深藍(lán)+暗像元，空間分辨率10 km)，為獲得較高時(shí)空分辨率的氣溶膠產(chǎn)品，將2種產(chǎn)品進(jìn)行融合，融合原則是以MOD04_3K數(shù)據(jù)為主，用MOD04_L2數(shù)據(jù)補(bǔ)充缺失區(qū)域。

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

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

(1)

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

2 研究方法

2.1 基本原理

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

(2)

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

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

y=0.49x+0.05

(3)

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

圖2 MODIS與WFI紅藍(lán)波段地表反射率散點(diǎn)圖

ρ3=0.500 4ρ1+0.046 8

(4)

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

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

圖3 MODIS與WFI藍(lán)波段地表反射率擬合圖

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

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

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

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

(5)

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

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

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

(6)

式中：ρ1為藍(lán)光波段值；ρ3為紅光波段值；ρ4為近紅外波段值。

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

3 結(jié)果與分析

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

3.1 氣溶膠反演方法對(duì)比與分析

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

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

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

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

圖6 WFI反演AOD與MODIS AOD空間分布對(duì)比圖

3.3 基于AERONET站點(diǎn)監(jiān)測(cè)數(shù)據(jù)的精度驗(yàn)證

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

EE=±(0.05+0.2·τα)

(7)

式中：τα為AERONET站點(diǎn)實(shí)測(cè)AOD值。

如圖7所示，本研究共得220個(gè)有效數(shù)據(jù)對(duì)，反演結(jié)果與AERONET站點(diǎn)實(shí)測(cè)值相關(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散點(diǎn)對(duì)比圖

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

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

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

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