朱 世 鑫，潘 竟 虎

(西北師范大學(xué)地理與環(huán)境科學(xué)學(xué)院，甘肅 蘭州 730070)

0 引言

大氣氣溶膠是懸浮在空氣中的分散顆粒，包括初級(jí)氣溶膠(以顆粒的形式直接排放到大氣中)和次級(jí)氣溶膠(由大氣中的主要污染物轉(zhuǎn)化而來(lái))，直徑為0.001～100 μm[1]。大氣氣溶膠的濃度變化直接影響人體健康與空氣質(zhì)量，并通過(guò)輻射強(qiáng)迫的直接或間接效應(yīng)影響地氣收支平衡和氣候變化[2，3]，因此，對(duì)氣溶膠的空間分布和變化進(jìn)行監(jiān)測(cè)非常重要。氣溶膠光學(xué)厚度(Aerosol Optical Depth,AOD)是描述氣溶膠特性的重要參數(shù)之一，可表示大氣的渾濁度[4]。傳統(tǒng)氣溶膠監(jiān)測(cè)多為地基監(jiān)測(cè)，然而由于站點(diǎn)布設(shè)有限，難以獲取大范圍的氣溶膠空間分布信息。近年來(lái)，利用具有高時(shí)效性和大尺度觀測(cè)能力的衛(wèi)星遙感技術(shù)可獲取空間分布連續(xù)的氣溶膠數(shù)據(jù)，在宏觀環(huán)境和污染分布監(jiān)測(cè)上具有較大潛力[5]。

AOD反演的關(guān)鍵是確定地表反射率信息，以便于實(shí)現(xiàn)地氣解耦[6]，反演算法包括暗目標(biāo)算法[7]、深藍(lán)算法[8,9]、結(jié)構(gòu)函數(shù)法[10]、偏振算法[11]、多角度算法[12]等。其中，結(jié)構(gòu)函數(shù)法需找到“清潔日”影像作為基準(zhǔn)且對(duì)幾何校正要求較高，偏振算法僅能應(yīng)用于反演細(xì)粒子氣溶膠，多角度算法需要特定的傳感器支持，3種算法的業(yè)務(wù)化應(yīng)用較為困難。近年來(lái)暗目標(biāo)算法不斷得到改進(jìn)與應(yīng)用[13,14]，但其仍具有明顯局限性，不僅需要短波紅外波段的支持，且不適用于亮像元地區(qū)(城市、沙漠等區(qū)域)。深藍(lán)算法可以反演出亮像元地區(qū)的AOD，但其精度低于暗目標(biāo)算法，而且需要外部的地表反射率產(chǎn)品，目前常使用MOD09A1數(shù)據(jù)，空間分辨率為500 m；由于單波段的反演需要具體的地表反射率柵格值才能實(shí)現(xiàn)地氣解耦，在反演過(guò)程中，如果表觀反射率數(shù)據(jù)的空間分辨率小于500 m，通常會(huì)將其重采樣為500 m[15-17]，最終得到500 m分辨率的反演結(jié)果，因此，深藍(lán)算法反演結(jié)果受外部地表反射率數(shù)據(jù)空間分辨率限制。如何能在不損失原有數(shù)據(jù)空間分辨率的前提下，同時(shí)反演出暗、亮像元地區(qū)氣溶膠光學(xué)厚度，是當(dāng)前氣溶膠光學(xué)厚度反演研究的一個(gè)重點(diǎn)。

Sentinel-3A衛(wèi)星搭載的OLCI(Ocean Land Colour Instrument)傳感器有21個(gè)通道，波段范圍為400～1 020 nm，輻射分辨率和信噪比較高，AOD反演潛力很好，且其空間分辨率為300 m，高于常使用的MODIS傳感器，但目前對(duì)于Sentinel-3A OLCI傳感器的氣溶膠反演研究還較少。Mei等[18]利用XBAER算法進(jìn)行OLCI數(shù)據(jù)的AOD反演，但其構(gòu)建地表反射率數(shù)據(jù)庫(kù)的過(guò)程限制了算法應(yīng)用。王峰等[19]采用紅藍(lán)通道比值算法，根據(jù)NASA發(fā)布的MODIS紅藍(lán)通道地表反射率關(guān)系式，經(jīng)過(guò)光譜轉(zhuǎn)換后獲取適用于Sentinel-3A OLCI傳感器的固定關(guān)系式，對(duì)中國(guó)臺(tái)灣暗像元地區(qū)(濃密植被區(qū))進(jìn)行AOD反演的結(jié)果較好，但該反演結(jié)果不包括亮像元地區(qū)，且根據(jù)NASA得到的固定關(guān)系式?jīng)]有考慮季節(jié)性和區(qū)域性差異，最終影響了反演精度。

本文基于現(xiàn)有理論，參考深藍(lán)算法的地表反射率庫(kù)獲取方法，并根據(jù)紅藍(lán)通道比值算法的思想，以中原城市群作為研究區(qū)，擬合不同NDVI等級(jí)下的地表反射率關(guān)系式，建立逐月的地表反射率關(guān)系庫(kù)，利用擬合得到的系數(shù)和截距反演AOD，以期消除外部地表反射率數(shù)據(jù)對(duì)反演結(jié)果空間分辨率的限制，獲取同時(shí)包括暗、亮像元的高空間分辨率AOD反演結(jié)果，并利用SONET(Sun-Sky Radiometer Observation Network)地基觀測(cè)數(shù)據(jù)檢驗(yàn)AOD反演結(jié)果，用MODIS AOD產(chǎn)品與反演結(jié)果進(jìn)行空間分布對(duì)比分析，以此檢驗(yàn)反演結(jié)果的精度和可靠性。

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

1.1 研究區(qū)概況

根據(jù)國(guó)家發(fā)改委2016年12月發(fā)布的《中原城市群發(fā)展規(guī)劃》，中原城市群是中國(guó)東西向經(jīng)濟(jì)發(fā)展的過(guò)渡地帶，涵蓋河南省18個(gè)地級(jí)市以及河北、山西、安徽、山東等省份的部分地級(jí)市。該地區(qū)資源壓力大，產(chǎn)業(yè)結(jié)構(gòu)以第二產(chǎn)業(yè)為主，近年來(lái)發(fā)展速度較快，人們?cè)谌粘Ｉ詈兔芗陨a(chǎn)中向城市大氣中排放的氣溶膠粒子日益增多，大氣環(huán)境污染問(wèn)題日漸嚴(yán)峻，尤其在污染物不易擴(kuò)散的秋冬季更為嚴(yán)重[20，21]，將其作為研究區(qū)具有代表性。研究區(qū)內(nèi)可用的SONET地基觀測(cè)站點(diǎn)有焦作站和嵩山站(圖1)。

圖1 中原城市群范圍和SONET 站點(diǎn)位置Fig.1 Scope and SONET site location of Central Plains Urban Agglomeration

1.2 數(shù)據(jù)與預(yù)處理

1.2.1 AOD反演所需數(shù)據(jù)

(1)MOD09A1數(shù)據(jù)。該數(shù)據(jù)是經(jīng)大氣校正后的地表反射率數(shù)據(jù)，選取8天內(nèi)質(zhì)量最佳的數(shù)據(jù)組合，包含MODIS傳感器的Band1-Band7波段，用于構(gòu)建地表反射率關(guān)系庫(kù)。由于OLCI傳感器與MODIS傳感器之間存在光譜響應(yīng)差異，因此需要利用實(shí)測(cè)光譜數(shù)據(jù)對(duì)MOD09A1數(shù)據(jù)進(jìn)行修正，而光譜庫(kù)在一定程度上可以解決缺少實(shí)測(cè)數(shù)據(jù)的問(wèn)題[22]。將ENVI光譜庫(kù)中的不同地物光譜重采樣至1 nm間隔，與MODIS和OLCI的紅、藍(lán)波段的光譜響應(yīng)函數(shù)分別進(jìn)行卷積計(jì)算，得到不同地物在MODIS或OLCI的紅、藍(lán)波段下的真實(shí)反射率(式(1))；對(duì)不同地物真實(shí)反射率進(jìn)行散點(diǎn)擬合，由式(2)模擬計(jì)算MODIS紅(藍(lán))波段與OLCI紅(藍(lán))波段的光譜轉(zhuǎn)換統(tǒng)計(jì)模型，獲取Sentinel-3A的地表反射率數(shù)據(jù)。本文采用最小值合成法處理每月的地表反射率數(shù)據(jù)，獲得每月地表反射率最小值，以去除云、陰影等影響，并將結(jié)果作為當(dāng)月的地表反射率真實(shí)值。

(1)

式中：Γ(λi)表示波長(zhǎng)為λi時(shí)的傳感器響應(yīng)函數(shù)；R(λi)為λi處的地物光譜值；Δλ采用0.0001微分值。

ROLCI=ARMODIS+B

(2)

式中：ROLCI和RMODIS分別表示由式(1)計(jì)算得到的OLCI和MODIS的真實(shí)地表反射率；A、B為對(duì)應(yīng)的擬合參數(shù)。

(2)Sentinel-3A數(shù)據(jù)。該數(shù)據(jù)來(lái)源于歐空局(https://scihub.copernicus.eu)，本文選取研究區(qū)晴空像元相對(duì)較多的L1B級(jí)影像數(shù)據(jù)進(jìn)行AOD反演，時(shí)間范圍為2019年3月至2020年2月；除各波段的輻射亮度數(shù)據(jù)外，輔助數(shù)據(jù)包括太陽(yáng)(衛(wèi)星)天頂角、太陽(yáng)(衛(wèi)星)方位角數(shù)據(jù)。利用SNAP軟件對(duì)原始輻射亮度數(shù)據(jù)和幾何數(shù)據(jù)進(jìn)行預(yù)處理，得到投影轉(zhuǎn)換后的表觀反射率數(shù)據(jù)和角度數(shù)據(jù)，并使用編程語(yǔ)言進(jìn)行數(shù)據(jù)裁剪及波段合成。

(3)DEM數(shù)據(jù)。該數(shù)據(jù)來(lái)源于地理空間數(shù)據(jù)云(http://www.gscloud.cn)的SRTMDEMUTM 90M數(shù)據(jù)集，用于對(duì)AOD反演結(jié)果進(jìn)行高程校正。為便于與Sentinel-3A數(shù)據(jù)匹配，將90 m的DEM數(shù)據(jù)重采樣為300 m。

1.2.2 驗(yàn)證數(shù)據(jù)

(1)SONET觀測(cè)數(shù)據(jù)。SONET是由中國(guó)科學(xué)院聯(lián)合中國(guó)多所高等院校和科研院所等單位在中國(guó)典型區(qū)域建立的觀測(cè)網(wǎng)，相比于集中分布于京津冀、臺(tái)灣等地區(qū)的AERONET(Aerosol Robotic Network)觀測(cè)數(shù)據(jù)，SONET站點(diǎn)在中國(guó)境內(nèi)分布更均勻，覆蓋全國(guó)多個(gè)省域，目前共有23個(gè)站點(diǎn)。此外，SONET具有與AERONET相當(dāng)?shù)臄?shù)據(jù)采集能力，且不確定性小于AERONET[23]。因此，本文選取研究區(qū)內(nèi)的焦作站和嵩山站的Level 1.5級(jí)觀測(cè)數(shù)據(jù)(http://www.sonet.ac.cn)，用于驗(yàn)證AOD反演結(jié)果的精度。由于反演結(jié)果是550 nm處的AOD，但站點(diǎn)觀測(cè)數(shù)據(jù)沒(méi)有550 nm處的AOD值，不便于直接對(duì)比，所以需要使用已有數(shù)據(jù)計(jì)算出550 nm處的AOD值作為驗(yàn)證數(shù)據(jù)。相較于Angstrom波長(zhǎng)指數(shù)法，通過(guò)二次多項(xiàng)式法計(jì)算得到的AOD插值數(shù)據(jù)結(jié)果更可靠[24]。該方法使用3個(gè)已知通道的AOD擬合未知參數(shù)，然后計(jì)算出550 nm處的AOD值(式(3))，本文選取675 nm、500 nm、440 nm 3個(gè)通道進(jìn)行計(jì)算。

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

(3)

式中：τλ為λnm處AOD值；a0、a1、a2為對(duì)應(yīng)的擬合參數(shù)。

(2)MOD04_3K數(shù)據(jù)。本文選用MOD04_3K AOD作為反演結(jié)果空間分布的檢驗(yàn)數(shù)據(jù)，其與Sentinel-3A衛(wèi)星過(guò)境時(shí)間相近，是NASA發(fā)布的Level 2級(jí)第6版氣溶膠數(shù)據(jù)，空間分辨率3 km，在濃密植被區(qū)反演精度較高。

2 研究方法

2.1 基本原理

假設(shè)大氣分子的分布與變化均一，且地表是均勻朗伯面，則衛(wèi)星接收到表觀反射率(ρTOA)的計(jì)算公式為：

ρTOA(us,uv,φ)=

(4)

式中：us、uv分別為太陽(yáng)天頂角和衛(wèi)星天頂角的余弦值；φ為相對(duì)方位角；Tg為總氣體分子透過(guò)率；ρ0為大氣路徑程輻射項(xiàng)；T(us)、T(uv)分別為大氣上行、下行透過(guò)率；ρs為地表反射率；S為大氣下界的半球反射率。將式(4)變換可得地表反射率的求解方程：

ρs(us,uv,φ)=

(5)

由式(4)、式(5)可知，表觀反射率中包含大氣和地表的信息，因此，AOD反演的前提是實(shí)現(xiàn)地氣解耦，去除地表信息，然后從大氣信息中得到AOD的反演結(jié)果。具體流程(圖2)為：利用輻射傳輸模型獲取預(yù)設(shè)氣溶膠類型條件下的參數(shù)查找表，結(jié)合表觀反射率與角度文件，迭代輸入不同AOD值得到相應(yīng)的大氣參數(shù)，并根據(jù)式(5)計(jì)算出對(duì)應(yīng)的紅藍(lán)波段地表反射率；然后與已知的真實(shí)地表反射率關(guān)系式進(jìn)行對(duì)比，當(dāng)二者最接近時(shí)，當(dāng)前AOD迭代值即為AOD反演值。

圖2 研究流程Fig.2 Flow chart of the study

2.2 查找表構(gòu)建

計(jì)算地表反射率所用的參數(shù)一般通過(guò)大氣輻射傳輸模型獲取，目前常使用6SV(Second Simulation of the Satellite Signal in the Solar Spectrum vector code)模型。但6SV 2.1模型中并沒(méi)有OLCI傳感器的光譜響應(yīng)函數(shù)，需要利用FORTRAN語(yǔ)言對(duì)模型源碼進(jìn)行修改，添加OLCI傳感器的光譜響應(yīng)支持。在構(gòu)建查找表時(shí)使用添加的波段，可減少由光譜響應(yīng)函數(shù)不同造成的誤差[25]。

若直接利用大氣輻射傳輸模型計(jì)算，速度緩慢，為提高反演效率，通常采用查找表的方式進(jìn)行AOD反演。本文在查找表的設(shè)置中，將太陽(yáng)天頂角和衛(wèi)星方位角的步長(zhǎng)設(shè)置為6°，相對(duì)方位角的步長(zhǎng)設(shè)置為12°，550 nm處AOD值的步長(zhǎng)設(shè)置為0.05，高程設(shè)置為0，大氣類型設(shè)置為中緯度冬季和中緯度夏季；此外，相較于其他氣溶膠模式，大陸型氣溶膠模式反演結(jié)果精度較高[26]，因此，本文將氣溶膠模式設(shè)置為大陸型氣溶膠。最終得到550 nm處AOD值和3種幾何角度在不同組合條件下的大氣參數(shù)查找表，參數(shù)包括：總氣體分子透過(guò)率(Tg)、大氣路徑程輻射項(xiàng)(ρ0)、總大氣透過(guò)率(T(us)×T(uv))、大氣下界的半球反射率(S)，對(duì)于查找表中不存在的幾何角度所對(duì)應(yīng)參數(shù)，通過(guò)線性插值的方式獲取。

2.3 云水掩膜

云像元的存在會(huì)對(duì)大氣參數(shù)的精度造成影響，因此云掩膜不充分會(huì)導(dǎo)致AOD反演結(jié)果出現(xiàn)較大誤差。常用的云掩膜方法需要借助熱紅外波段，從云層與地表存在溫差的角度剔除云像元。由于OLCI傳感器沒(méi)有熱紅外波段，同時(shí)考慮到根據(jù)表觀反射率及其空間變化特征設(shè)置合適的閾值，可有效地實(shí)現(xiàn)云檢測(cè)[27]，本文采用如下針對(duì)OLCI傳感器的快速云檢測(cè)方法(圖3)：在厚云區(qū)域，OA3通道中的表觀反射率ρTOA3較高，利用此性質(zhì)，以0.35作為閾值識(shí)別厚云像元[19]；根據(jù)薄云區(qū)域的非均一性，當(dāng)ρTOA3的3×3像元標(biāo)準(zhǔn)差大于0.002時(shí)，或當(dāng)ρTOA6的3×3像元標(biāo)準(zhǔn)差大于0.05時(shí)，判別為云像元；OA13通道中的表觀反射率ρTOA13隨云頂高度的增加而增加，而OA12通道中的表觀反射率ρTOA12對(duì)云頂高度的依賴性非常弱[28]，因此二者在云像元的反射率比值大于非云像元的反射率比值，設(shè)定比值大于0.315時(shí)為云像元；此外，對(duì)初步的云識(shí)別結(jié)果邊界進(jìn)行5×5像元擴(kuò)充處理，以去除云邊緣和云陰影的不確定性。從圖4可以看出，由于缺少熱紅外波段的支持，卷積云區(qū)域出現(xiàn)輕微漏檢，但云掩膜算法整體運(yùn)行穩(wěn)定，能較準(zhǔn)確地識(shí)別出云像元區(qū)域。

注:STD3×3表示3×3像元標(biāo)準(zhǔn)差，ρTOAi表示OAi通道中的表觀反射率。

利用OA8和OA17通道的表觀反射率計(jì)算NDVI，將NDVI小于0的區(qū)域判定為水體像元進(jìn)行水掩膜。

2.4 地表反射率關(guān)系庫(kù)構(gòu)建

本文基于紅藍(lán)通道比值算法的思想，通過(guò)擬合不同NDVI下的紅藍(lán)波段地表反射率線性關(guān)系，建立地表反射率關(guān)系庫(kù)，利用固定關(guān)系式，根據(jù)系數(shù)和截距進(jìn)行反演，避免了計(jì)算的地表反射率與真實(shí)地表反射率之間的比較，從而實(shí)現(xiàn)高空間分辨率的氣溶膠光學(xué)厚度反演。

首先，將OLCI表觀反射率重采樣到500 m，以便與地表反射率數(shù)據(jù)匹配。考慮到濃密植被區(qū)與非濃密植被區(qū)的差異，通過(guò)表觀反射率計(jì)算NDVI，根據(jù)NDVI劃分等級(jí)分別進(jìn)行擬合：濃密植被區(qū)的紅藍(lán)波段線性關(guān)系較穩(wěn)定，參考王峰等[19]的研究，將NDVI大于0.45的像元判斷為濃密植被區(qū)；非濃密植被區(qū)由于植被稀疏程度不同，地表反射率關(guān)系較為復(fù)雜，因此，每0.05的NDVI設(shè)置一個(gè)級(jí)別，擬合不同NDVI等級(jí)下的紅藍(lán)波段地表反射率關(guān)系。最后，將擬合出的地表反射率關(guān)系式應(yīng)用到原始的Sentinel-3A數(shù)據(jù)中，反演出300 m分辨率的AOD。擬合過(guò)程如下：假設(shè)短時(shí)間內(nèi)地表反射率不發(fā)生變化或變化較小，以最小值合成法合成的地表反射率數(shù)據(jù)作為當(dāng)月的真實(shí)地表反射率，使用每月的多景待反演影像計(jì)算NDVI，去除云、水像元后取NDVI均值，以保證有足夠的像元參與地表反射率關(guān)系的擬合；然后讀取不同NDVI下對(duì)應(yīng)位置的紅藍(lán)波段地表反射率像元值，以NDVI等級(jí)分類并分別進(jìn)行擬合，得到一階線性方程式，建立地表反射率關(guān)系庫(kù)?？紤]到擬合關(guān)系的季節(jié)性差異，該地表反射率關(guān)系庫(kù)是逐月建立的，表1為2020年2月紅藍(lán)波段地表反射率關(guān)系擬合結(jié)果。

圖4 云識(shí)別效果與影像對(duì)比Fig.4 Comparison between cloud recognition effect and remote sensing image

表1 2020年2月地表反射率關(guān)系庫(kù)Table 1 Relationship database of land surface reflectivity in February 2020

2.5 AOD反演

利用經(jīng)過(guò)云、水掩膜處理后的300 m表觀反射率數(shù)據(jù)計(jì)算出NDVI，根據(jù)地表反射率關(guān)系庫(kù)確定該NDVI下使用的擬合關(guān)系式；然后以0.001作為初始AOD進(jìn)行計(jì)算，步長(zhǎng)β設(shè)置為0.05，根據(jù)預(yù)處理后的表觀反射率和幾何數(shù)據(jù)，代入式(5)迭代計(jì)算不同AOD條件下的紅、藍(lán)波段地表反射率ρ8和ρ4。由于表觀反射率會(huì)隨著AOD的增加而增加，且波長(zhǎng)更短的波段對(duì)氣溶膠增加更敏感[29]，因此，隨著AOD的增加，(ρ4-b)/ρ8的結(jié)果會(huì)不斷減小，當(dāng)(ρ4-b)/ρ8小于系數(shù)a時(shí)停止迭代，并根據(jù)迭代前后的數(shù)據(jù)插值獲取初步反演結(jié)果(圖5)。

圖5 AOD反演算法流程Fig.5 Flow chart of AOD retrieval algorithm

2.6 高程校正

單個(gè)查找表無(wú)法兼顧研究區(qū)內(nèi)所有高程信息，建立多個(gè)高程條件下的查找表會(huì)影響計(jì)算效率，因此，本文以高程為0設(shè)置查找表，對(duì)AOD初步反演結(jié)果進(jìn)行高程校正[30]。在目標(biāo)高度為Z(單位為km，標(biāo)準(zhǔn)大氣高度為8.5 km)時(shí)，大氣分子的光學(xué)厚度可表示為：

τλ(z=Z)=τλ(z=0)exp(-Z/8.5)

(6)

2.7 精度驗(yàn)證

本文使用SONET站點(diǎn)數(shù)據(jù)對(duì)反演所得OLCI AOD進(jìn)行精度驗(yàn)證，評(píng)價(jià)指標(biāo)包括決定系數(shù)(R2)、均方根誤差(RMSE)、期望誤差(EE)[31]，計(jì)算公式如下：

(7)

EE=±(0.05+0.2×τ0)

(8)

式中：τ為反演AOD值；τ0為站點(diǎn)AOD值；n為驗(yàn)證樣本數(shù)。

3 結(jié)果與驗(yàn)證

3.1 空間分布一致性檢驗(yàn)

利用上述算法獲得研究區(qū)內(nèi)300 m空間分辨率的OLCI AOD數(shù)據(jù)，并與MOD04_3K AOD數(shù)據(jù)相比較，選擇的影像類型包括：1)研究區(qū)AOD值區(qū)間較大，出現(xiàn)區(qū)域性霧霾天氣(圖6a和圖6d)。OLCI反演結(jié)果與MOD04_3K AOD均可有效識(shí)別出研究區(qū)中部的焦作市、新鄉(xiāng)市和鄭州市與北部的邢臺(tái)市、聊城市、邯鄲市和安陽(yáng)市的區(qū)域性霧霾天氣，但OLCI AOD的空間分辨率更高，可清晰地識(shí)別出中部地區(qū)霧霾天氣的細(xì)節(jié)信息；在城市群北部的邢臺(tái)市，本文估算結(jié)果的有效像元比MODIS產(chǎn)品更多，MODIS產(chǎn)品在邢臺(tái)市AOD高值區(qū)出現(xiàn)很多缺失像元，可能是由于重霧霾天氣下的大氣散射作用導(dǎo)致該區(qū)域反射率偏大，被誤判為亮像元，因此MODIS在反演AOD時(shí)剔除了這些像元。2)研究區(qū)AOD值區(qū)間較小，空氣質(zhì)量較好(圖6b和圖6e)。OLCI反演結(jié)果與MOD04_3K AOD仍保持一致，AOD值較低的區(qū)域主要分布于西部秦嶺等地區(qū)和西北部太行山區(qū)，在東部平原地區(qū)能明顯看出城市區(qū)域AOD值較高,表明人類活動(dòng)對(duì)氣溶膠分布具有較顯著影響。3)MOD04_3K AOD產(chǎn)品在晴空像元區(qū)域出現(xiàn)大量數(shù)據(jù)缺失(圖6c和圖6f)。北方地區(qū)冬季地表植被覆蓋相對(duì)較少、地表反射率高，因此基于暗目標(biāo)算法的MOD04_3K AOD產(chǎn)品在研究區(qū)內(nèi)數(shù)據(jù)缺失較多，而本文反演的OLCI AOD仍可獲取大范圍的AOD分布，從而反映出研究區(qū)的霧霾分布特征。

圖6 OLCI AOD與MOD04 AOD空間分布對(duì)比Fig.6 Comparison of spatial distribution between OLCI AOD and MOD04 AOD

綜上可知，OLCI AOD與MOD04_3K AOD之間高值、低值區(qū)域相互對(duì)應(yīng)，具有較好的空間分布一致性，但本文基于地表反射率關(guān)系庫(kù)反演的AOD的空間分辨率更高，展現(xiàn)的細(xì)節(jié)信息更多，并且空間覆蓋有效范圍也更大。

3.2 SONET站點(diǎn)驗(yàn)證

首先，為使驗(yàn)證結(jié)果更合理，對(duì)過(guò)境時(shí)間前后的兩個(gè)觀測(cè)結(jié)果進(jìn)行線性插值，以此作為驗(yàn)證數(shù)據(jù)；然后進(jìn)行驗(yàn)證數(shù)據(jù)與反演結(jié)果的相關(guān)分析，并計(jì)算出精度驗(yàn)證評(píng)價(jià)指標(biāo)，評(píng)估本文反演結(jié)果的精度(圖7)。從圖7可以看出，本文反演結(jié)果與SONET觀測(cè)結(jié)果之間相關(guān)性較好，決定系數(shù)R2達(dá)到0.9099，但部分反演結(jié)果存在高估現(xiàn)象，RMSE低至0.1909，偏差整體較小，有67.44%的OLCI AOD在期望誤差(EE)區(qū)間內(nèi)，進(jìn)一步說(shuō)明本研究的反演結(jié)果整體反演精度和可靠性較高。

圖7 OLCI AOD與SONET AOD精度驗(yàn)證Fig.7 Precision verification of OLCI AOD based on SONET AOD

根據(jù)圖8a可以發(fā)現(xiàn)嵩山站點(diǎn)的反演結(jié)果較好，OLCI AOD與SONET的觀測(cè)值和變化趨勢(shì)較為接近，但當(dāng)AOD接近0.1時(shí)，反演結(jié)果出現(xiàn)明顯低估，這可能是由查找表的構(gòu)建和插值方法所致[22]。焦作站點(diǎn)的反演結(jié)果(圖8b)與觀測(cè)值相比整體偏高，尤其是3月的反演結(jié)果高估現(xiàn)象較為明顯,主要原因是該月份與焦作站點(diǎn)同等級(jí)NDVI的像元數(shù)量較少，影響了關(guān)系庫(kù)的擬合精度，進(jìn)而影響到反演結(jié)果，此外，地表反射率的數(shù)據(jù)質(zhì)量也會(huì)影響反演結(jié)果。

圖8 OLCI AOD與SONET AOD時(shí)序?qū)Ρ菷ig.8 Comparison of time series between OLCI AOD and SONET AOD

3.3 誤差分析

3.3.1 擬合誤差影響 本文AOD反演是基于紅藍(lán)通道比值算法中固定系數(shù)與截距的思想。系數(shù)和截距是影響反演結(jié)果的重要因素，有必要進(jìn)一步分析其差異性導(dǎo)致的誤差，因此本文基于6SV模型定量模擬二者對(duì)AOD反演造成的誤差。紅藍(lán)波段地表反射率比值多處于1.65～2.25之間[32]，因此定量模擬中假設(shè)比值為2.0，探究不同反射率地表下系數(shù)和截距所造成的誤差，參數(shù)設(shè)定與分析結(jié)果如圖9所示，其中，OZA，SZA，RAA，ρ8分別表示衛(wèi)星天頂角、太陽(yáng)天頂角、相對(duì)方位角和紅光波段地表反射率，大氣類型為中緯度夏季。

圖9 擬合誤差影響分析Fig.9 Analysis of the influence of fitting error

由圖9a可知，0.01的截距誤差大約會(huì)造成0.1的反演誤差，不同地表反射率下的誤差較接近。由圖9b可知，隨著地表反射率增大，系數(shù)誤差引起的反演誤差也增大，且偏差較大：當(dāng)ρ8小于0.15時(shí)，0.1的系數(shù)誤差約造成0.1的反演誤差；當(dāng)ρ8大于0.15時(shí)，0.1的系數(shù)誤差至少造成0.2的反演誤差。因此，隨著地表反射率增大，擬合誤差對(duì)反演結(jié)果的影響會(huì)更大。

3.3.2 誤差來(lái)源分析 擬合誤差是反演結(jié)果最重要的影響因素，此外，以下因素也在一定程度上影響反演結(jié)果。1)數(shù)據(jù)質(zhì)量與數(shù)量：在使用地表反射率數(shù)據(jù)計(jì)算地表反射率關(guān)系式時(shí)，采用最小值合成法去除云、陰影的影響，但該方法不一定能完全剔除受影響像元，從而導(dǎo)致擬合誤差，進(jìn)而影響反演結(jié)果；此外，部分月份內(nèi)某些NDVI等級(jí)所對(duì)應(yīng)的像元數(shù)量不足，也會(huì)導(dǎo)致擬合質(zhì)量不佳，從而造成反演誤差。2)氣溶膠類型：本文選擇了普適性更強(qiáng)的大陸型氣溶膠，未考慮氣溶膠組分所引起的變化。隨著AOD值的增加，氣溶膠類型的影響會(huì)越來(lái)越重要[33]。3)地表反射率變化：本文假設(shè)地表反射率在短期內(nèi)無(wú)變化或變化較小，然而在自然或人為影響下(如冰雪消融、農(nóng)作物收割等),地表反射率會(huì)發(fā)生較大改變，造成反演誤差[2]。4)太陽(yáng)位置和傳感器位置的差異：本文雖然已對(duì)MODIS地表反射率進(jìn)行光譜轉(zhuǎn)換以消除傳感器不同產(chǎn)生的誤差，但太陽(yáng)位置、傳感器位置等差異也會(huì)對(duì)反演結(jié)果產(chǎn)生一定影響。

4 結(jié)論

本文考慮區(qū)域性、下墊面不同及季節(jié)性變化，假設(shè)短期內(nèi)地表反射率不變或變化較小，通過(guò)逐月擬合不同下墊面條件下的紅藍(lán)波段地表反射率關(guān)系式，建立地表反射率關(guān)系庫(kù)，實(shí)現(xiàn)了高空間分辨率的AOD反演，得到以下結(jié)論：1)紅藍(lán)通道比值算法依據(jù)固定關(guān)系式進(jìn)行反演，建立關(guān)系庫(kù)可以避免與具體地表反射率柵格值的對(duì)比，不受外部地表反射率數(shù)據(jù)的限制，從而獲取更高空間分辨率的反演結(jié)果。2)與SONET站點(diǎn)數(shù)據(jù)的驗(yàn)證結(jié)果顯示，本文算法不僅能夠獲取更高空間分辨率的反演結(jié)果，而且反演質(zhì)量也較高，決定系數(shù)R2達(dá)到0.9099，RMSE低至0.1909，67.44%的OLCI AOD落在預(yù)期誤差區(qū)間(EE)內(nèi)。3)系數(shù)和截距是影響反演結(jié)果的重要因素，0.1的系數(shù)誤差會(huì)造成至少0.1的反演誤差，0.01的截距誤差大約會(huì)造成0.1的反演誤差。

本文借助真實(shí)地表反射率數(shù)據(jù)建立了逐月的地表反射率關(guān)系庫(kù)，由于小區(qū)域范圍內(nèi)并不能保證每一個(gè)NDVI等級(jí)內(nèi)都有足夠的像元進(jìn)行關(guān)系擬合，從而導(dǎo)致部分月份反演誤差相對(duì)較大，未來(lái)可以考慮引入歷史同期數(shù)據(jù)，保證關(guān)系庫(kù)建立所需的數(shù)據(jù)數(shù)量及質(zhì)量達(dá)到要求；此外，還需進(jìn)一步研究太陽(yáng)位置、傳感器位置等不同導(dǎo)致的反演結(jié)果誤差。