基于地表反射率關系庫的Sentinel-3A AOD反演
——以中原城市群為例

朱 世 鑫，潘 竟 虎

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

0 引言

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

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

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

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

1 研究區和數據

1.1 研究區概況

根據國家發改委2016年12月發布的《中原城市群發展規劃》，中原城市群是中國東西向經濟發展的過渡地帶，涵蓋河南省18個地級市以及河北、山西、安徽、山東等省份的部分地級市。該地區資源壓力大，產業結構以第二產業為主，近年來發展速度較快，人們在日常生活和密集性生產中向城市大氣中排放的氣溶膠粒子日益增多，大氣環境污染問題日漸嚴峻，尤其在污染物不易擴散的秋冬季更為嚴重[20，21]，將其作為研究區具有代表性。研究區內可用的SONET地基觀測站點有焦作站和嵩山站(圖1)。

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

1.2 數據與預處理

1.2.1 AOD反演所需數據

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

(1)

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

ROLCI=ARMODIS+B

(2)

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

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

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

1.2.2 驗證數據

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

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

(3)

式中：τλ為λnm處AOD值；a0、a1、a2為對應的擬合參數。

(2)MOD04_3K數據。本文選用MOD04_3K AOD作為反演結果空間分布的檢驗數據，其與Sentinel-3A衛星過境時間相近，是NASA發布的Level 2級第6版氣溶膠數據，空間分辨率3 km，在濃密植被區反演精度較高。

2 研究方法

2.1 基本原理

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

ρTOA(us,uv,φ)=

(4)

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

ρs(us,uv,φ)=

(5)

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

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

2.2 查找表構建

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

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

2.3 云水掩膜

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

注:STD3×3表示3×3像元標準差，ρTOAi表示OAi通道中的表觀反射率。

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

2.4 地表反射率關系庫構建

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

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

圖4 云識別效果與影像對比Fig.4 Comparison between cloud recognition effect and remote sensing image

表1 2020年2月地表反射率關系庫Table 1 Relationship database of land surface reflectivity in February 2020

2.5 AOD反演

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

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

2.6 高程校正

單個查找表無法兼顧研究區內所有高程信息，建立多個高程條件下的查找表會影響計算效率，因此，本文以高程為0設置查找表，對AOD初步反演結果進行高程校正[30]。在目標高度為Z(單位為km，標準大氣高度為8.5 km)時，大氣分子的光學厚度可表示為：

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

(6)

2.7 精度驗證

本文使用SONET站點數據對反演所得OLCI AOD進行精度驗證，評價指標包括決定系數(R2)、均方根誤差(RMSE)、期望誤差(EE)[31]，計算公式如下：

(7)

EE=±(0.05+0.2×τ0)

(8)

式中：τ為反演AOD值；τ0為站點AOD值；n為驗證樣本數。

3 結果與驗證

3.1 空間分布一致性檢驗

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

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

綜上可知，OLCI AOD與MOD04_3K AOD之間高值、低值區域相互對應，具有較好的空間分布一致性，但本文基于地表反射率關系庫反演的AOD的空間分辨率更高，展現的細節信息更多，并且空間覆蓋有效范圍也更大。

3.2 SONET站點驗證

首先，為使驗證結果更合理，對過境時間前后的兩個觀測結果進行線性插值，以此作為驗證數據；然后進行驗證數據與反演結果的相關分析，并計算出精度驗證評價指標，評估本文反演結果的精度(圖7)。從圖7可以看出，本文反演結果與SONET觀測結果之間相關性較好，決定系數R2達到0.9099，但部分反演結果存在高估現象，RMSE低至0.1909，偏差整體較小，有67.44%的OLCI AOD在期望誤差(EE)區間內，進一步說明本研究的反演結果整體反演精度和可靠性較高。

圖7 OLCI AOD與SONET AOD精度驗證Fig.7 Precision verification of OLCI AOD based on SONET AOD

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

圖8 OLCI AOD與SONET AOD時序對比Fig.8 Comparison of time series between OLCI AOD and SONET AOD

3.3 誤差分析

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

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

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

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

4 結論

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

本文借助真實地表反射率數據建立了逐月的地表反射率關系庫，由于小區域范圍內并不能保證每一個NDVI等級內都有足夠的像元進行關系擬合，從而導致部分月份反演誤差相對較大，未來可以考慮引入歷史同期數據，保證關系庫建立所需的數據數量及質量達到要求；此外，還需進一步研究太陽位置、傳感器位置等不同導致的反演結果誤差。