牛曉君，唐家奎?，張自力，崔林麗，幸煒鵬，宋 藝

(1 中國(guó)科學(xué)院大學(xué)資源與環(huán)境學(xué)院， 北京 100049； 2 浙江省環(huán)境監(jiān)測(cè)中心，杭州 310015； 3 上海市氣象科學(xué)研究所，上海 200030)

近年來(lái)，隨著經(jīng)濟(jì)和城市化迅猛發(fā)展，中國(guó)霧霾污染日趨嚴(yán)重,霧霾遙感監(jiān)測(cè)有迫切的應(yīng)用需求。傳統(tǒng)的霧霾監(jiān)測(cè)手段主要是利用稀疏的地基站點(diǎn)進(jìn)行PM2.5測(cè)量，僅能獲得觀測(cè)點(diǎn)處的空氣污染狀況，對(duì)于大范圍的霧霾監(jiān)測(cè)，衛(wèi)星遙感具有較強(qiáng)優(yōu)勢(shì)[1]。中國(guó)霧霾災(zāi)害往往具有季節(jié)性、區(qū)域性特點(diǎn)，具有較多可見光譜段的極軌衛(wèi)星已經(jīng)初步展示霧霾監(jiān)測(cè)的能力，以TERRA和AQUA衛(wèi)星為例，每天能實(shí)現(xiàn)1～2次全球覆蓋觀測(cè)，搭載的MODIS傳感器獲取的氣溶膠光學(xué)厚度(aerosol optical depth, AOD)產(chǎn)品是目前常用霧霾監(jiān)測(cè)的指標(biāo)[2]。然而,霧霾的主要成分氣溶膠的生命周期短、時(shí)空變化強(qiáng)，每日1～2次的極軌衛(wèi)星觀測(cè)無(wú)法滿足對(duì)其性質(zhì)、動(dòng)態(tài)變化等特性觀測(cè)的需求，尤其對(duì)于霧霾的發(fā)生、發(fā)展、擴(kuò)散遷移、消散等過(guò)程的監(jiān)測(cè)，這需要更高頻次的觀測(cè)。靜止衛(wèi)星具有高時(shí)間分辨率的特點(diǎn)，可以對(duì)地面進(jìn)行定點(diǎn)凝視觀測(cè)，觀測(cè)頻率可達(dá)小時(shí)級(jí)甚至分鐘級(jí)，因此，靜止衛(wèi)星對(duì)于霧霾動(dòng)態(tài)過(guò)程的監(jiān)測(cè)具有明顯優(yōu)勢(shì)[3]。

目前，霧霾(PM2.5)衛(wèi)星遙感監(jiān)測(cè)工作主要是基于搭載衛(wèi)星上的氣溶膠監(jiān)測(cè)傳感器開展的，在衛(wèi)星遙感反演氣溶膠的基礎(chǔ)上，再進(jìn)一步進(jìn)行霧霾監(jiān)測(cè)分析[1]。因此，氣溶膠光學(xué)厚度反演成為霧霾遙感監(jiān)測(cè)的關(guān)鍵環(huán)節(jié)。國(guó)內(nèi)外已有學(xué)者開展了利用靜止衛(wèi)星監(jiān)測(cè)氣溶膠光學(xué)厚度及霧霾災(zāi)害的研究工作。毛節(jié)泰和劉莉[4]開展用GMS5衛(wèi)星資料反演湖面上空氣溶膠的可行性研究，結(jié)果表明，通過(guò)選取合適的湖面反射率,衛(wèi)星反演的氣溶膠光學(xué)厚度和地面光度計(jì)遙感的月均值相對(duì)誤差不超過(guò) 30%。張軍華等[5]利用GMS衛(wèi)星遙感中國(guó)大陸25個(gè)湖面上空氣溶膠光學(xué)厚度，結(jié)果指出,利用GMS衛(wèi)星反演湖面光學(xué)厚度可以提供大陸上空氣溶膠的信息。高玲等[6]利用日本靜止衛(wèi)星MTSAT數(shù)據(jù)采用最小衛(wèi)星反射率反演地表反射率的方法反演氣溶膠光學(xué)厚度，并與AERONET觀測(cè)及MODIS氣溶膠光學(xué)厚度產(chǎn)品進(jìn)行比較，結(jié)果表明，MTSAT反演的氣溶膠光學(xué)厚度產(chǎn)品可以反映大氣氣溶膠光學(xué)厚度的日變化信息。任通等[3]探討利用中國(guó)風(fēng)云2C靜止衛(wèi)星可見光資料反演氣溶膠光學(xué)厚度的數(shù)值方法，研究結(jié)果表明，在東亞地區(qū)利用風(fēng)云2C可見光資料反演的AOD產(chǎn)品可以反映氣溶膠的分布，但算法存在不同地區(qū)AOD值的高估和低估現(xiàn)象。陳健等[7]應(yīng)用韓國(guó)靜止衛(wèi)星GOCI數(shù)據(jù)，采用國(guó)際深藍(lán)算法進(jìn)行長(zhǎng)江三角洲地區(qū)氣溶膠光學(xué)厚度的遙感反演研究，發(fā)現(xiàn)基于靜止衛(wèi)星GOCI數(shù)據(jù)的氣溶膠產(chǎn)品具有監(jiān)測(cè)氣溶膠日變化的能力。Moulin 等[8]利用歐空局靜止衛(wèi)星Meteosat進(jìn)行氣溶膠反演，并計(jì)算1983—1994年期間由撒哈拉沙漠輸送到大西洋和地中海上空的沙塵量。Wang等[9]利用靜止衛(wèi)星GOES8反演大西洋區(qū)域的沙塵氣溶膠，并與地面觀測(cè)試驗(yàn)進(jìn)行比對(duì)，結(jié)果表明，高時(shí)間分辨率的靜止衛(wèi)星可以有效彌補(bǔ)極軌衛(wèi)星全球氣溶膠遙感的不足。Prados等[10]評(píng)估NOAA/NESDIS利用靜止衛(wèi)星GEOS可見光通道反演的氣溶膠/煙塵產(chǎn)品(GASP)，通過(guò)與地基AERONET和MODIS官方產(chǎn)品比對(duì)，表明反演結(jié)果在北美地區(qū)具有較好的精度。

日本氣象廳于2014年7月發(fā)射新一代靜止氣象衛(wèi)星Himawari-8，Himawari-8可實(shí)現(xiàn)時(shí)間分辨率為10 min/次的高頻次對(duì)地觀測(cè)，引起廣泛關(guān)注。目前已有利用該衛(wèi)星傳感器應(yīng)用于天氣預(yù)報(bào)、海面溫度、灰霾檢測(cè)等方面的相關(guān)研究[11-13]。Himawari-8搭載的AHI傳感器具有對(duì)氣溶膠敏感的可見光通道。目前已有學(xué)者利用Himawari-8數(shù)據(jù)進(jìn)行氣溶膠遙感和霧霾監(jiān)測(cè)的嘗試：Daisaku等[14]利用Himawari-8雙通道(可見光0.64 μm和近紅外0.86 μm波段)數(shù)據(jù)進(jìn)行氣溶膠反演試驗(yàn)；Sekiyama等[15]將Himawari-8反演的氣溶膠光學(xué)厚度用于全球氣溶膠預(yù)測(cè)模型Kalman濾波系統(tǒng)，并預(yù)測(cè)了東亞的一次沙塵暴。葛邦宇等[16]利用Himawari-8的可見光和近紅外波段進(jìn)行了氣溶膠反演，并驗(yàn)證了用暗像元法反演葵花8號(hào)衛(wèi)星數(shù)據(jù)得到的結(jié)果精度及方法的可行性。Wang等[17]利用葵花衛(wèi)星官方氣溶膠產(chǎn)品，采用線性混合效應(yīng)模型對(duì)2015年7月至2017年5月京津冀地區(qū)的PM2.5相關(guān)分析建模，結(jié)果表明葵花衛(wèi)星能夠有效用于PM2.5監(jiān)測(cè)。

本文提出一種利用葵花衛(wèi)星反演氣溶膠光學(xué)厚度的算法，并將其結(jié)果與AERONET、MODIS等觀測(cè)的氣溶膠光學(xué)厚度值進(jìn)行對(duì)比，對(duì)反演結(jié)果進(jìn)行誤差分析。最后以氣溶膠光學(xué)厚度為霧霾指標(biāo)，對(duì)北京區(qū)域一次霧霾過(guò)程進(jìn)行反演分析。

1 衛(wèi)星及數(shù)據(jù)

葵花8號(hào)衛(wèi)星(Himawari 8)是由日本氣象廳于2014年發(fā)射的氣象衛(wèi)星(http:∥www.eorc.jaxa.jp/ptree/userguide.html)，經(jīng)過(guò)性能檢測(cè)后，于2015年正式投入使用?？?號(hào)衛(wèi)星的星下點(diǎn)位于東經(jīng)140.7°的赤道上空，距地高度約為35 800 km，搭載有AHI成像儀，該傳感器共有3個(gè)可見光通道和13個(gè)近紅外及紅外通道，最高分辨率可達(dá)0.5 km。AHI全盤掃描周期為10 min，如選定時(shí)間對(duì)特定區(qū)域進(jìn)行掃描，可達(dá)2.5 min觀測(cè)周期。相對(duì)于極軌衛(wèi)星，葵花8號(hào)靜止衛(wèi)星具有空間分辨率高、觀測(cè)頻率高、觀測(cè)范圍廣的應(yīng)用優(yōu)勢(shì)，目前葵花衛(wèi)星數(shù)據(jù)已經(jīng)應(yīng)用于天氣、海洋、灰霾、火點(diǎn)監(jiān)測(cè)等。由于AHI傳感器設(shè)置了對(duì)氣溶膠較為敏感的可見光波段，因此在氣溶膠反演方面潛力巨大，已有葵花衛(wèi)星官方產(chǎn)品包括AOD(aerosol optical depth)、AMV(atmospheric motion vector)、CSR(clear sky radiance)、HCAI(high-resolution cloud analysis information)。本文將采用葵花衛(wèi)星AHI傳感器可見光紅、藍(lán)波段數(shù)據(jù)(分辨率分別為0.5和1 km)進(jìn)行氣溶膠光學(xué)厚度反演，并結(jié)合AERONET站點(diǎn)觀測(cè)數(shù)據(jù)、MODIS氣溶膠產(chǎn)品以及葵花衛(wèi)星官方氣溶膠產(chǎn)品(空間分辨率為5 km)進(jìn)行算法精度、空間分布和空間覆蓋率方面的比對(duì)驗(yàn)證。

本研究采用的地面監(jiān)測(cè)驗(yàn)證數(shù)據(jù)來(lái)自于國(guó)際AERONET觀測(cè)站。AERONET (AErosol RObotic NETwork)是一個(gè)國(guó)際性的氣溶膠地基觀測(cè)網(wǎng)站(https:∥aeronet.gsfc.nasa.gov)，由NASA和PHOTONS在全球范圍共建。主要提供3個(gè)質(zhì)量等級(jí)的數(shù)據(jù)：level 1.0為未進(jìn)行云檢測(cè)的數(shù)據(jù)，level 1.5為進(jìn)行了云去除的數(shù)據(jù)，level 2.0為進(jìn)行了云去除和質(zhì)量檢測(cè)的數(shù)據(jù)。本文選用的對(duì)比站點(diǎn)為北京地區(qū)的Beijing、Beijing-CAMS、Beijing-PKU這3個(gè)站點(diǎn)，位置分布見圖1。由于2017年尚未發(fā)布level 2.0數(shù)據(jù)，驗(yàn)證數(shù)據(jù)所以采用level 1.5作為驗(yàn)證數(shù)據(jù)。

(底圖是2017年5月17上午10時(shí)葵花衛(wèi)星數(shù)據(jù)真彩色合成圖，分辨率為1 km)圖1 北京地區(qū)AERONET站點(diǎn)Beijing、Beijing-CAMS、Beijing-PKU分布圖Fig.1 Location map of AERONET sites of Beijing， Beijing-CAMS, and Beijing-PKU in Beijing area

本文同時(shí)采用目前流行的極軌衛(wèi)星MODIS官方氣溶膠產(chǎn)品進(jìn)行比對(duì)。MODIS是搭載在Terra和Aqua兩顆衛(wèi)星上的中分辨率成像光譜儀，其數(shù)據(jù)可以從官網(wǎng)(https:∥ladsweb.modaps.eosdis.nasa.gov)免費(fèi)獲得。本文采用MOD04暗像元法和深藍(lán)算法結(jié)合的氣溶膠產(chǎn)品作為交叉驗(yàn)證數(shù)據(jù)，其空間分辨率為10 km。

2 葵花衛(wèi)星反演AOD算法

郭強(qiáng)等[18]提出利用MODIS可見光波段反演陸地氣溶膠光學(xué)厚度的算法，該算法無(wú)需事先獲得地表反射率先驗(yàn)知識(shí)或者如經(jīng)典暗像元法通過(guò)近紅外通道與可見光通道之間比值關(guān)系估算地表反射率，而是利用在較短時(shí)間內(nèi)地物可見光波段反射率比值變化較小的特性進(jìn)行反演。本文將該算法改進(jìn)應(yīng)用于葵花衛(wèi)星數(shù)據(jù)，算法原理如下：

根據(jù)Vermote等[19]研究結(jié)果，在地表朗伯體、大氣水平均一的假設(shè)條件下，衛(wèi)星平臺(tái)觀測(cè)到的大氣頂反射率(表觀反射率)可以表示為

ρTOA(μs,μs,φ)=

(1)

式中：ρ0是大氣路徑輻射項(xiàng)等效反射率，ρs是地表反射率，T(μs)、T(μv)分別是入射方向和觀測(cè)方向的大氣透過(guò)率，S為大氣下界面的半球反射率。在公式中T(μs)、T(μv)總是同時(shí)出現(xiàn)，因此將T(μs)T(μv)作為一個(gè)參數(shù)考慮。未知數(shù)ρ0,S,T(μs)T(μv)與AOD、氣溶膠粒子單次散射率、相函數(shù)、不對(duì)稱指數(shù)等有關(guān)。在實(shí)際反演過(guò)程中，首先針對(duì)試驗(yàn)區(qū)域確定一個(gè)氣溶膠模型(因本文研究區(qū)為北京，反演中默認(rèn)選擇大陸型氣溶膠模式)，然后用6S、MODTRAN等輻射傳輸模型模擬在不同觀測(cè)幾何情況下AOD與上述3個(gè)參數(shù)之間的對(duì)應(yīng)關(guān)系，建立查找表。若已知地表反射率和表觀反射率，便可通過(guò)查找表計(jì)算得到AOD。求解式(1)的困難之處就在于很難獲得準(zhǔn)確的地表反射率值。

根據(jù)已有的研究，較短時(shí)間內(nèi)(排除降雨、降雪等改變地表反射率情況)，同一地物的不同觀測(cè)幾何條件，可見光波段地表反射率比值變化較小，可忽略不計(jì)[20]，因此，針對(duì)同一地物兩次觀測(cè)，可以近似表達(dá)為

(2)

式中：r為地表反射率，m、n分別表示兩個(gè)不同的觀測(cè)幾何條件，b1、b2表示不同的波段。

根據(jù)式(2)可以導(dǎo)出

(3)

通過(guò)式(3)可以得出，對(duì)于同一地物在兩個(gè)觀測(cè)幾何條件下，分別對(duì)應(yīng)的兩個(gè)波段反射率的比值相同。經(jīng)典暗像元法采用較低反照率地表(暗像元)存在可見光波段(0.47、0.66 μm)與短波近紅外波段(2.13 μm)反射率比值經(jīng)驗(yàn)關(guān)系(即r0.47 μm=0.25r2.13 μm、r0.66 μm=0.5r2.13 μm)。因此，暗像元法的波段比值關(guān)系是存在于假設(shè)地物為朗伯體及地物光譜反射特性不隨時(shí)間變化的情況下，是本文式(3)的特例。

因此，對(duì)于葵花8 號(hào)衛(wèi)星AHI傳感器，紅、藍(lán)波段地表反射率之間可得如下關(guān)系

(4)

(5)

目前葵花衛(wèi)星尚無(wú)可靠的官方地表反射率產(chǎn)品發(fā)布，無(wú)法直接得到地表反射率比值k。MODIS官方地表反射率產(chǎn)品(MOD09)是目前國(guó)際上較為流行認(rèn)可采用的產(chǎn)品，已有學(xué)者采用MOD09用于估算其他衛(wèi)星的地表反射率[21-24]。

本文參考已有研究成果[23-24]，采用奇異值分解(SVD)方法實(shí)現(xiàn)應(yīng)用MODIS產(chǎn)品(MOD09)估算葵花衛(wèi)星數(shù)據(jù)地表反射率，進(jìn)而計(jì)算獲得葵花衛(wèi)星可見光波段地表反射率比值。Sayer等[23]利用ASTER、USGS地物光譜庫(kù)，分別模擬不同地物在MODIS以及AATSR觀測(cè)時(shí)的光譜反射率，利用基于奇異值分解(SVD)的方法，模擬建立MODIS和AATSR兩種傳感器的不同光譜響應(yīng)之間的轉(zhuǎn)換關(guān)系，從而實(shí)現(xiàn)利用MODIS的地表反射率模擬計(jì)算AATSR的地表反射率。郭強(qiáng)[24]同樣采用SVD方法應(yīng)用MODIS地表反射率產(chǎn)品計(jì)算HJ-1 A/1B CCD可見光波段地表反射率產(chǎn)品，進(jìn)而用于氣溶膠反演。前人研究結(jié)果表明，SVD分解方法可以很好地減少不同傳感器之間因?yàn)椴ǘ雾憫?yīng)函數(shù)的差異而產(chǎn)生的地表反射率誤差，能夠用于不同傳感器觀測(cè)的地表反射率之間的轉(zhuǎn)換。

奇異值分解(singular value decomposition, SVD)是一種統(tǒng)計(jì)的方法，可以用來(lái)提取一組數(shù)據(jù)的特征。假設(shè)存在一個(gè)m×n階矩陣M，那么矩陣M可以表示為

M=UΣV*.

(6)

本文利用Herold等[25]建立的光譜數(shù)據(jù)庫(kù)、MODIS光譜響應(yīng)函數(shù)、Himawri-8光譜響應(yīng)函數(shù)三者共同構(gòu)建矩陣M，選取其中光譜曲線125條，包括綠色植被、枯萎植被、裸土、水體以及屋頂、瀝青道路、水泥道路、停車場(chǎng)、鐵軌等共26種典型地物，用以代表城市及周邊大部分區(qū)域的典型地物類型。矩陣每行代表一個(gè)波段(前3行為MODIS、后3行為Himawari-8)，每一列代表一種光譜，矩陣M為6×125階矩陣。對(duì)矩陣M進(jìn)行奇異值分解后得到矩陣U，從矩陣U中選取前n個(gè)奇異向量構(gòu)成矩陣S：

S=[U1,U2,…，Un]

(7)

取矩陣S的前3行構(gòu)成S的子矩陣SM，后3行構(gòu)成子矩陣SH。對(duì)于矩陣M，第i個(gè)列向量的前3個(gè)元素代表第i種光譜的MODIS 反射率，記作向量Vim。同理，后三個(gè)元素代表H8反射率，記作向量Vih，對(duì)于向量Vim，構(gòu)建一個(gè)列向量Ci，存在

(8)

求解后可以得到模擬MODIS地表反射率Vim與模擬葵花衛(wèi)星反射率Vih的轉(zhuǎn)換關(guān)系

(9)

本文采用上述SVD方法，利用MOD09地表反射率產(chǎn)品估算H8可見光波段的地表反射率產(chǎn)品，進(jìn)而計(jì)算獲得紅藍(lán)波段比值k。實(shí)際反演求解式(5)時(shí)，采用基于6S輻射傳輸模型的價(jià)值函數(shù)迭代求解方法，即在給定氣溶膠模式、紅藍(lán)波段比值k情況下，對(duì)地表反射率和氣溶膠光學(xué)厚度分別以一定步長(zhǎng)間隔取值，代入6S輻射傳輸模型進(jìn)行計(jì)算，分別得到模擬紅、藍(lán)波段衛(wèi)星表觀反射率，建立模擬的衛(wèi)星表觀反射率和觀測(cè)衛(wèi)星表觀反射率之間的絕對(duì)損失函數(shù)(loss function)如下

(10)

3 AOD反演驗(yàn)證

3.1 AOD反演的精度驗(yàn)證

為了驗(yàn)證本文提出算法的反演精度，以北京作為研究區(qū)域，選取2017年5月中旬的一個(gè)霧霾過(guò)程開展研究。首先從日本官方網(wǎng)站下載到葵花衛(wèi)星L1級(jí)可見光波段數(shù)據(jù)，并進(jìn)行輻射校正、幾何投影和云檢測(cè)等預(yù)處理工作，由于云檢測(cè)不是本文研究重點(diǎn)，本文采取目視解譯方法嚴(yán)格剔除云像元。

葵花衛(wèi)星可以獲得每個(gè)10 min的高頻監(jiān)測(cè)數(shù)據(jù)，為驗(yàn)證算法精度，本文首先針對(duì)3個(gè)AERONET站點(diǎn)位置上空給出每10 min的葵花衛(wèi)星反演結(jié)果，并分別與AERONET觀測(cè)、葵花衛(wèi)星官方產(chǎn)品(L2，空間分辨率為5 km)、MODIS官方產(chǎn)品(空間分辨率為10 km)進(jìn)行對(duì)比驗(yàn)證?？臻g匹配時(shí)，本文算法反演結(jié)果和MODIS官方產(chǎn)品、葵花8號(hào)官方產(chǎn)品均采用對(duì)應(yīng)AERONET站點(diǎn)所在經(jīng)緯度位置為中心的3×3網(wǎng)格平均值；時(shí)間匹配時(shí)，選取靠近衛(wèi)星觀測(cè)最近時(shí)刻的AERONET站點(diǎn)觀測(cè)數(shù)據(jù)用于對(duì)比。結(jié)果對(duì)比如圖2。

圖2 AOD反演結(jié)果(H8retrieval_AOD)與AERONET(Beijing_AOD、 Beijing_CAMS_AOD、 Beijing_PKU_AOD代表站點(diǎn)觀測(cè)值)、葵花官方產(chǎn)品(H8ARP_AOD)及MODIS官方產(chǎn)品(MOD/MYD04)的對(duì)比驗(yàn)證Fig.2 Comparison of AOD retrieval results(H8retrieval_AOD)with AERONET observation data(Beijing_AOD, Beijing_CAMS_AOD, and Beijing_PKU_AOD), Himawari-8 official products(H8ARP_AOD), and MODIS official products(MOD/MYD04)

從圖2可以看到，總體上，利用本文算法反演得到的AOD值與AERONET站點(diǎn)數(shù)據(jù)波動(dòng)趨勢(shì)較為一致，并且值的大小也較為接近，在16時(shí)之后，反演結(jié)果與觀測(cè)值相比誤差較大，且均大于觀測(cè)值。可能因?yàn)?6時(shí)之后太陽(yáng)天頂角和相對(duì)方位角均較大的緣故導(dǎo)致。

統(tǒng)計(jì)可得，與AERONET 3站點(diǎn)Beijing、Beijing_CAMS、Beijing_PKU觀測(cè)結(jié)果相比，本文反演結(jié)果的絕對(duì)誤差平均值分別為0.06、0.05、0.03，相對(duì)誤差平均值分別為12%、10%、7%?？ü俜疆a(chǎn)品僅在10時(shí)至15時(shí)之間有值，絕對(duì)誤差平均值分別為0.12、0.11、0.07，相對(duì)誤差平均值分別為20%、18%、14%，總體上看，葵花官方產(chǎn)品與AERONET站點(diǎn)數(shù)據(jù)相比在上午偏差較大，且存在低估現(xiàn)象。TERRA和AQUA雙星MODIS官方氣溶膠產(chǎn)品上下午各有一個(gè)觀測(cè)值，絕對(duì)誤差平均值分別為0.27、0.25、0.33，相對(duì)誤差平均值分別為48%、46%、67%，MODIS官方產(chǎn)品誤差最大，且存在明顯高估現(xiàn)象。比較而言，不難看出，本文提出的算法反演精度明顯優(yōu)于葵花官方產(chǎn)品和MODIS官方產(chǎn)品，且在可反演的日時(shí)間范圍上長(zhǎng)于葵花官方產(chǎn)品。

本文反演算法的誤差來(lái)源可能包括以下幾個(gè)方面；1)經(jīng)驗(yàn)設(shè)置的氣溶膠類型參數(shù)，本文采取的是大陸型氣溶膠模式，而對(duì)于北京區(qū)域并不一定符合實(shí)際狀況；2)地表反射率比值的估算精度，本文應(yīng)用MOD09地表反射率產(chǎn)品模擬計(jì)算獲得地表反射率比值，并非真實(shí)的葵花地表反射率產(chǎn)品計(jì)算的比值；3)目視解譯方法進(jìn)行的云檢測(cè)可能遺漏薄云像元；4)優(yōu)化求解、迭代求解時(shí)精度設(shè)置等。本文的反演案例中，AOD均小于1，對(duì)于更為嚴(yán)重的霧霾過(guò)程，如AOD?1時(shí)，需要進(jìn)一步研究，因6S模型在能見度低于5 km時(shí)的輻射傳輸計(jì)算可靠性降低，理論上分析，本文算法對(duì)于大氣能見度低于5 km時(shí)的AOD反演誤差會(huì)進(jìn)一步增大，因此較適用于能見度大于5 km時(shí)的霧霾狀況。

上述衛(wèi)星數(shù)據(jù)反演結(jié)果與地面觀測(cè)站點(diǎn)對(duì)比驗(yàn)證的誤差，除受反演結(jié)果本身精度影響外，還與對(duì)比驗(yàn)證時(shí)的時(shí)空匹配精度有關(guān)?？臻g上，相對(duì)于10 km分辨率的MODIS產(chǎn)品和5 km的葵花產(chǎn)品，本文反演的1 km空間分辨率產(chǎn)品更容易與地面站點(diǎn)匹配，另外對(duì)比時(shí)采用的是反演結(jié)果的3×3格網(wǎng)均值，因此在對(duì)比時(shí)空間匹配更具優(yōu)勢(shì)。時(shí)間上，相比于葵花衛(wèi)星10 min一次觀測(cè)頻率，MODIS只有上下午2次過(guò)境，與地面觀測(cè)的時(shí)間匹配更難吻合，時(shí)空匹配精度低將進(jìn)一步增大驗(yàn)證誤差。

3.2 AOD反演的空間分布驗(yàn)證

為進(jìn)一步分析驗(yàn)證本文提出算法反演的空間分布的準(zhǔn)確性和優(yōu)勢(shì)，本文將葵花衛(wèi)星的反演結(jié)果與NASA官網(wǎng)獲得的TERRA星上MODIS氣溶膠產(chǎn)品(MOD04)進(jìn)行對(duì)比，為了對(duì)比時(shí)相上的一致性，盡量選擇MODIS過(guò)境時(shí)刻接近的葵花衛(wèi)星反演結(jié)果。對(duì)比結(jié)果見圖3。

(a)5月17日上午10時(shí)55分MOD04氣溶膠產(chǎn)品圖像；(b)5月17日上午11時(shí)葵花衛(wèi)星反演結(jié)果圖像;(c)5月18日上午11時(shí)44分MOD04氣溶膠產(chǎn)品圖像;(d)5月18日上午11時(shí)50分葵花衛(wèi)星反演結(jié)果圖像圖3 葵花衛(wèi)星反演結(jié)果圖像與MOD04氣溶膠產(chǎn)品圖像對(duì)比Fig.3 Comparison between retrieval AOD image and MODIS official AOD products image

圖3中，顏色由藍(lán)到紅表示氣溶膠光學(xué)厚度從低到高。由兩顆衛(wèi)星的AOD反演結(jié)果對(duì)比可以看出，總體上氣溶膠光學(xué)厚度的空間分布較為一致，整體呈現(xiàn)東南高，西北低的狀況。在東南部城市區(qū)域，MOD04值明顯比本文算法的反演結(jié)果要高，這與前述圖2中MOD04明顯高于本文算法反演結(jié)果和AERONET站點(diǎn)觀測(cè)值的對(duì)比結(jié)果是一致的。另外，由于本文的反演結(jié)果空間分辨率(1 km)明顯優(yōu)于MODIS官方產(chǎn)品(10 km)，可以看出葵花衛(wèi)星反演結(jié)果較MODIS產(chǎn)品空間分布細(xì)節(jié)清晰。另外也可以看出，MOD04產(chǎn)品存在較多的無(wú)值空白區(qū)域，這是由于MODIS產(chǎn)品采用的暗像元法局限性導(dǎo)致，而葵花衛(wèi)星反演算法不受亮地表限制，反演結(jié)果有效區(qū)域明顯大于MOD04。因此，在時(shí)間分辨率、空間分辨率以及應(yīng)用范圍的角度，葵花衛(wèi)星反演結(jié)果均優(yōu)于MODIS官方應(yīng)用產(chǎn)品。

進(jìn)一步地，將本文算法反演結(jié)果與葵花衛(wèi)星官方產(chǎn)品進(jìn)行對(duì)比?？ü俜綒馊苣z產(chǎn)品反演也是采用暗像元法，產(chǎn)品目前仍處于實(shí)驗(yàn)階段。目前官方提供的兩種質(zhì)量級(jí)別的氣溶膠產(chǎn)品，即L2(10 min級(jí))和L3(1 h級(jí))兩級(jí)數(shù)據(jù)，空間分辨率均為5 km，其中L3級(jí)數(shù)據(jù)是L2級(jí)數(shù)據(jù)的小時(shí)平均值。由于算法受太陽(yáng)天頂角限制較大，在中國(guó)大陸地區(qū)每日10時(shí)后才有反演產(chǎn)品，另外除云覆蓋影響外，由于算法限制，葵花官方氣溶膠產(chǎn)品在中國(guó)覆蓋率極低，大部分時(shí)間只有中國(guó)東部沿海地區(qū)有數(shù)據(jù)，除云覆蓋影響外，也許與暗像元法可應(yīng)用的區(qū)域局限性有關(guān)。如圖4所示，兩個(gè)級(jí)別產(chǎn)品在中國(guó)均有較大范圍的空白無(wú)效區(qū)域，其中L3產(chǎn)品是L2數(shù)據(jù)經(jīng)過(guò)后續(xù)篩選處理獲得，精度質(zhì)量控制高于L2級(jí)數(shù)據(jù)，但時(shí)間分辨率由10 min級(jí)降為1 h級(jí)，有效監(jiān)測(cè)空間范圍進(jìn)一步縮小，無(wú)法滿足中國(guó)大范圍區(qū)域監(jiān)測(cè)應(yīng)用的需求。

圖4 2017年5月17日11時(shí)葵花衛(wèi)星官方氣溶膠產(chǎn)品(ARP)L2級(jí)(a)和L3級(jí)(b)Fig.4 Himawari-8 satellite official aerosol products (ARP) of level 2 (a) and level 3 (b) at 11 a.m. on May 17 2017

為了兼顧比較L2級(jí)(更大空間覆蓋范圍)和L3級(jí)(更好的精度質(zhì)量控制)數(shù)據(jù)的優(yōu)勢(shì)，本文同時(shí)選取葵花官方產(chǎn)品每日11時(shí)至14時(shí)整點(diǎn)時(shí)刻的L2和L3級(jí)數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證(見圖5)，反演結(jié)果圖像時(shí)間分別為：11、12、13、14整點(diǎn)時(shí)刻。

圖5 葵花衛(wèi)星反演結(jié)果圖像與葵花衛(wèi)星官方氣溶膠產(chǎn)品圖像對(duì)比Fig.5 Comparison between retrieval AOD image and Himawari-8 official AOD products image

總體上三者AOD的空間分布趨勢(shì)較為一致，2個(gè)日期均在北京東南部出現(xiàn)AOD高值區(qū)。比較而言，在2個(gè)日期里，L3級(jí)產(chǎn)品圖像均在西北部出現(xiàn)大范圍空白無(wú)效區(qū)域，這是因?yàn)長(zhǎng)3產(chǎn)品為L(zhǎng)2數(shù)據(jù)后續(xù)篩選處理的結(jié)果，說(shuō)明葵花官方產(chǎn)品在優(yōu)化產(chǎn)品精度時(shí)犧牲了產(chǎn)品空間覆蓋率；L2級(jí)產(chǎn)品圖像上零星分布無(wú)效空白像元，在整幅圖像上呈明顯斑點(diǎn)狀，但空間覆蓋范圍明顯優(yōu)于L3級(jí)數(shù)據(jù)。綜合空間覆蓋范圍和空間分布細(xì)節(jié)兩個(gè)方面，本文算法反演結(jié)果均為最優(yōu)?？梢钥闯觯谒惴▍^(qū)域適應(yīng)性方面，本文算法優(yōu)于葵花官方產(chǎn)品采用的暗像元法，同時(shí)，由于本文算法反演結(jié)果空間分辨率(1 km)明顯高于葵花官方產(chǎn)品(5 km)，能夠更好滿足較高空間分辨率AOD產(chǎn)品的需要。

4 霧霾過(guò)程監(jiān)測(cè)

霧霾的主要成因是對(duì)流層氣溶膠，常用表征指標(biāo)是PM2.5濃度或空氣質(zhì)量指數(shù)AQI，而大氣氣溶膠屬性常用表征指標(biāo)是氣溶膠光學(xué)厚度AOD，其定義為從大氣底面到大氣頂層的垂直氣柱內(nèi)氣溶膠消光系數(shù)的積分。由于PM2.5濃度或AQI常常采用的是近地面測(cè)量數(shù)據(jù)，而衛(wèi)星遙感AOD是整層大氣氣溶膠消光系數(shù)的積分，因此衛(wèi)星遙感監(jiān)測(cè)霧霾常常采用AOD作為間接大氣污染程度指標(biāo)或應(yīng)用AOD估算近地面PM2.5濃度或AQI，諸多研究表明氣溶膠光學(xué)厚度與大氣污染程度呈正相關(guān)性，可以定性反映地面大氣污染狀況[25-26]。如前所述，本文對(duì)北京區(qū)域2017年5月17日、18日進(jìn)行小時(shí)級(jí)氣溶膠AOD反演，反演結(jié)果如圖6所示。

圖6 2017年5月17日和18日8時(shí)—17時(shí)小時(shí)級(jí)時(shí)間序列氣溶膠AOD反演結(jié)果圖像Fig.6 Hourly time series AOD retrieval images from 8 a.m. to 17 p.m. on May 17 and 18, 2017

本文獲得35個(gè)北京市空氣質(zhì)量監(jiān)測(cè)站點(diǎn)2017年5月17日8時(shí)至17時(shí)及5月18日8時(shí)至16時(shí)(17時(shí)數(shù)據(jù)缺失)的PM2.5監(jiān)測(cè)數(shù)據(jù)。依據(jù)站點(diǎn)的經(jīng)緯度地理位置以及監(jiān)測(cè)數(shù)據(jù)采集時(shí)間，分別與本文的AOD反演結(jié)果圖像進(jìn)行時(shí)空匹配，提取出站點(diǎn)位置處的PM2.5與AOD的數(shù)據(jù)對(duì)，進(jìn)行相關(guān)性分析，如圖7所示。

本文反演的AOD與地面監(jiān)測(cè)的PM2.5之間直接對(duì)比的皮爾森相關(guān)系數(shù)為0.502，滿足統(tǒng)計(jì)學(xué)上99%置信度的要求，這說(shuō)明PM2.5與AOD之間存在一定程度的線性正相關(guān)，但是相關(guān)性不是很強(qiáng)。前人的研究結(jié)果表明如果進(jìn)一步對(duì)反演的AOD進(jìn)行濕度訂正與垂直訂正，有望進(jìn)一步提高AOD和PM2.5的相關(guān)性[26-27]。即便如此，亦可進(jìn)一步說(shuō)明氣溶膠光學(xué)厚度與大氣污染程度呈一定的正相關(guān)性，可以用來(lái)定性反映地面大氣污染狀況，作為城市地面顆粒物濃度監(jiān)測(cè)的一個(gè)有效補(bǔ)充手段。

為驗(yàn)證本文算法對(duì)霧霾過(guò)程監(jiān)測(cè)的可行性，以AOD為霧霾程度監(jiān)測(cè)指標(biāo)，分析北京區(qū)域此次霧霾過(guò)程的時(shí)空動(dòng)態(tài)變化。從這兩日反演AOD圖像中，可以清楚地看到，總體上北京區(qū)域的霧霾分布呈現(xiàn)東南污染重，西北污染相對(duì)輕的趨勢(shì)?？傮w上看，5月18日霧霾污染相比于5月17日更為嚴(yán)重。經(jīng)查詢中國(guó)空氣質(zhì)量在線監(jiān)測(cè)分析平臺(tái)，歷史數(shù)據(jù)記錄顯示5月17日北京區(qū)域?yàn)橹卸任廴?，空氣質(zhì)量指數(shù)AQI日均值為184，5月18日為重度污染，AQI日均值為201。說(shuō)明本文反演的AOD的日變化趨勢(shì)與官方大氣環(huán)境監(jiān)測(cè)數(shù)據(jù)結(jié)果基本一致。

在時(shí)間序列上，本文采用的葵花衛(wèi)星數(shù)據(jù)反演結(jié)果具有明顯的時(shí)間分辨率優(yōu)勢(shì)，從5月17日8時(shí)至5月18日17時(shí)，每隔一小時(shí)的反演結(jié)果顯示AOD空間分布一直在變化，尤其18日相對(duì)17日而言，氣溶膠空間分布變化更快，說(shuō)明霧霾時(shí)空分布的高度動(dòng)態(tài)變化性特征。時(shí)間序列上，可以看出，5月17日8時(shí)起北京城區(qū)開始輕度霧霾，到下午17時(shí)，城中心霧霾加重，且在城南鄰近河北區(qū)域出現(xiàn)霧霾嚴(yán)重區(qū)域。5月18日上午8時(shí)起北京城區(qū)域再次出現(xiàn)輕度霧霾，且從上午11時(shí)起，霧霾開始明顯加重，至下午13點(diǎn)達(dá)到霧霾最重，然后至下午17時(shí)，霧霾開始減弱。另外也可以看出，隨著時(shí)間變化，霧霾污染中心具有逐漸向北京城東北方向擴(kuò)散移動(dòng)的趨勢(shì)。

從這兩日反演結(jié)果時(shí)間序列圖像中，可以清楚地看到霧霾污染中心的移動(dòng)方向以及加重或減輕的動(dòng)態(tài)變化狀況，證明霧霾動(dòng)態(tài)變化強(qiáng)，具有輸送、擴(kuò)散等特點(diǎn)。在實(shí)際應(yīng)用中，使用傳統(tǒng)的極軌衛(wèi)星數(shù)據(jù)如MODIS等獲得的氣溶膠數(shù)據(jù)已無(wú)法滿足對(duì)霧霾動(dòng)態(tài)高頻次監(jiān)測(cè)的要求。目前霧霾監(jiān)測(cè)預(yù)報(bào)主要使用地面監(jiān)測(cè)站點(diǎn)，空間上分布稀疏，難以進(jìn)行大范圍監(jiān)測(cè)和預(yù)報(bào)。霧霾的高空間分辨率、高時(shí)間分辨率以及大范圍同時(shí)監(jiān)測(cè)的需求，使得靜止氣象衛(wèi)星有望成為未來(lái)霧霾監(jiān)測(cè)較好的選擇。本文研究表明，利用葵花8號(hào)衛(wèi)星高時(shí)空分辨率可以進(jìn)行每1 h甚至每10 min級(jí)的反演，能夠動(dòng)態(tài)描繪出污染物的擴(kuò)散運(yùn)動(dòng)變化軌跡，如能結(jié)合當(dāng)時(shí)的氣象狀況，有望成為研究污染物的溯源和擴(kuò)散路徑追蹤的有效手段。

5 結(jié)論

本文提出一種利用葵花8號(hào)衛(wèi)星可見光波段反演高時(shí)空分辨率氣溶膠光學(xué)厚度的算法，并通過(guò)與國(guó)際AERONET 3個(gè)站點(diǎn)Beijing、Beijing-CAMS、Beijing-PKU站觀測(cè)數(shù)據(jù)，葵花官方氣溶膠產(chǎn)品和MODIS官方氣溶膠產(chǎn)品對(duì)比分析，驗(yàn)證算法的可行性。結(jié)果表明，本文提出算法的反演結(jié)果，無(wú)論在反演精度、有效反演空間范圍、有效反演時(shí)間范圍、空間分辨率等方面均為最優(yōu)。該算法不受地表反射類型的限制，可應(yīng)用于亮地表地區(qū)，突破了暗像元法的局限性。相比于MODIS氣溶膠數(shù)據(jù)和葵花官方氣溶膠數(shù)據(jù)，該算法的反演結(jié)果具有高覆蓋度，高精度、高頻次、高分辨率的特點(diǎn)，應(yīng)用本文算法可以獲得每10 min級(jí)、空間分辨率為1 km的氣溶膠光學(xué)厚度日動(dòng)態(tài)變化反演結(jié)果，更有利于對(duì)霧霾動(dòng)態(tài)變化的監(jiān)測(cè)。最后通過(guò)對(duì)北京區(qū)域2017年5月17日、18日小時(shí)級(jí)氣溶膠AOD反演結(jié)果的時(shí)空分析，探討北京區(qū)域此次霧霾過(guò)程的時(shí)空動(dòng)態(tài)變化，從這兩日反演結(jié)果時(shí)間序列圖像中，可以清楚地看到霧霾污染中心的移動(dòng)方向以及加重或減輕的動(dòng)態(tài)變化狀況，說(shuō)明本文提出的應(yīng)用葵花衛(wèi)星監(jiān)測(cè)霧霾過(guò)程的方法具有較好應(yīng)用潛力。存在的問(wèn)題及未來(lái)工作包括：1)云像元剔除由現(xiàn)在的目視解譯改成利用自動(dòng)云檢測(cè)算法實(shí)現(xiàn)動(dòng)態(tài)云掩膜；2)目前僅依據(jù)經(jīng)驗(yàn)采用大陸型氣溶膠模式，未來(lái)可根據(jù)實(shí)際情況采用自定義氣溶膠模式，進(jìn)一步提高AOD反演精度；3)本文的反演案例中，AOD均小于1，對(duì)于更為嚴(yán)重的霧霾過(guò)程，如AOD?1時(shí)，需要未來(lái)進(jìn)一步研究，理論上分析，本文算法較適用于能見度大于5 km時(shí)的霧霾狀況；4)目前計(jì)算地表反射率比值方法是基于MOD09產(chǎn)品的奇異值分解(SVD)模擬計(jì)算，未來(lái)將探索比較其他地表反射率比值的計(jì)算方法，進(jìn)一步提高精度；5)針對(duì)霧霾的監(jiān)測(cè)，本文直接應(yīng)用葵花衛(wèi)星反演的AOD對(duì)霧霾進(jìn)行定性分析，未做近地面AOD訂正處理，未來(lái)將開展反演AOD的濕度訂正和垂直訂正等后處理，獲得近地面AOD后進(jìn)一步估算顆粒物PM2.5、PM10濃度等研究。

感謝NASA官方網(wǎng)站提供的MODIS數(shù)據(jù)產(chǎn)品、日本氣象廳官方網(wǎng)站提供的葵花8號(hào)衛(wèi)星數(shù)據(jù)及國(guó)際AERONET站北京Beijing、Beijing-CAMS、Beijing-PKU 3站的數(shù)據(jù)支持。