田信鵬，高志強(qiáng)，劉強(qiáng)，王德，王躍啟

1.中國科學(xué)院煙臺海岸帶研究所 海岸帶環(huán)境過程與生態(tài)修復(fù)重點(diǎn)實(shí)驗(yàn)室,煙臺 264003;

2.中國科學(xué)院煙臺海岸帶研究所 山東省海岸帶環(huán)境過程重點(diǎn)實(shí)驗(yàn)室,煙臺 264003;

3.鵬成國家實(shí)驗(yàn)室,深圳 518055

1 引 言

氣溶膠光學(xué)厚度AOD（Aerosol Optical Depth）作為描述氣溶膠光學(xué)特性的重要參數(shù)之一，準(zhǔn)確獲取其時(shí)空分布信息是當(dāng)前大氣科學(xué)的研究熱點(diǎn)。衛(wèi)星遙感反演是獲取區(qū)域/全球尺度氣溶膠信息的有效技術(shù)手段，而地表反射率估算作為影響反演精度的主要因素之一（楊磊庫等，2022），降低地表參數(shù)的不確定性正是算法優(yōu)化的重要方面。經(jīng)過幾十年的發(fā)展，依據(jù)反演過程中采取的地表反射率估算策略差異產(chǎn)生了不同算法，如波段經(jīng)驗(yàn)比值（Kaufman 等，1997；Levy 等，2013）、多天反射率數(shù)據(jù)合成（Hsu 等，2013；田信鵬等，2018）、多源傳感器聯(lián)合反演（Tang 等，2005；Vermote 等，2007）、機(jī)器學(xué)習(xí)算法（賈臣 等，2020；Su 等，2020）等，氣溶膠遙感已形成了一個(gè)豐富的研究體系。

利用濃密植被地表反射率在紅、藍(lán)波段與短波紅外具有較穩(wěn)定的比值關(guān)系，且短波紅外受大氣影響相對較小的特點(diǎn)，Kaufman 等（1997）提出了適用于濃密植被地區(qū)的DT（Dark Target）算法。Levy 等（2013）基于DT 算法，結(jié)合散射角及短波紅外歸一化植被指數(shù)等，實(shí)現(xiàn)了MODIS氣溶膠產(chǎn)品（MOD04）業(yè)務(wù)化生產(chǎn)，目前已發(fā)展到C6.1 版本。葛邦宇等（2018）將DT 算法應(yīng)用于Himawari?8 衛(wèi)星，實(shí)現(xiàn)了京津冀地區(qū)AOD 反演。DT 算法能夠以較高精度獲取濃密植被地表上空AOD 信息，為提高氣溶膠衛(wèi)星遙感的空間覆蓋度，國內(nèi)外學(xué)者提出了多種適用于非濃密植被地表的反演算法。Hsu等（2004，2013）提出了適用于高亮地區(qū)的DB（Deep Blue） 算 法；孫 林 等（2016） 利 用 多 年MODIS 地表反射率產(chǎn)品構(gòu)建地表反射率先驗(yàn)數(shù)據(jù)集，提出了HARLS（High?resolution Aerosol Retrieval algorithm based ona prioriLand Surface reflectance database support）算法，并將其用于GF?1 PMS 等傳感器；Lyapustin 等（2018）提出的MAIAC（Multiangle implementation of atmospheric correction）算法是近年來發(fā)展起來的另一較成熟的反演算法，該算法采用動態(tài)時(shí)間序列數(shù)據(jù)實(shí)現(xiàn)地氣分離，能同時(shí)獲取濃密植被及高反射地表上空氣溶膠信息，基于該算法實(shí)現(xiàn)了MODIS 1 km 分辨率氣溶膠產(chǎn)品（MCD19A2）的生產(chǎn)。

以上DT、DB 及HARLS 算法假設(shè)地表為各向同性反射的朗伯體，在反演過程中使用的地表反射率僅與時(shí)間和位置有關(guān)。這種假設(shè)在草原等區(qū)域影響不大，但在方向反射特征明顯、地表結(jié)構(gòu)復(fù)雜的城市地區(qū)，會存在地表參數(shù)刻畫不確定性大的問題（Wei 等，2019），致使AOD 反演精度較低。由于不同地表類型，其反射率變化呈現(xiàn)不同的時(shí)空變異性，且地表方向反射特性使地表信息確定變得更加復(fù)雜和困難（Gatebe 和King，2016）。近年來，國內(nèi)外諸多學(xué)者在地表反射估算方面做了許多努力，Diner 等（2005）利用MISR傳感器多角度多通道觀測的特點(diǎn)，用經(jīng)驗(yàn)正交函數(shù)描述地表信息，實(shí)現(xiàn)了考慮地表方向反射的AOD 反演；Shi 等（2017）基于比值梯度優(yōu)化算法，實(shí)現(xiàn)了AATSR（Advanced Along Track Scanning Radiometer）傳感器AOD 及地表反射率同步反演。地表方向反射可用雙向反射分布函數(shù)BRDF（Bidirectional Reflectance Distribution Function） 刻畫，Xue等（2009）通過MODIS BRDF/Albedo 產(chǎn)品（MCD43）增加約束條件，基于AATSR傳感器星下點(diǎn)及前向兩個(gè)角度觀測數(shù)據(jù)，提出了協(xié)同反演算法；Yang 等（2014）利用非朗伯輻射傳輸前向模型，同時(shí)反演獲取了AOD 和?ngstrom 指數(shù)信息；Tian 等（2018） 結(jié)合半經(jīng)驗(yàn)核驅(qū)動模型估算地表反射率，實(shí)現(xiàn)了城市地區(qū)AOD 反演；She 等（2019）利用MCD43 提供的核權(quán)重信息，實(shí)現(xiàn)了Himawari?8 AHI 傳感器的AOD 和地表反射率聯(lián)合反演。

由于城市下墊面空間異質(zhì)性的影響，其地表反射特征刻畫困難，導(dǎo)致AOD 反演結(jié)果存在空間連續(xù)性差、分辨率低及精度具有較大不確定性等問題。Tian 等（2021a）通過耦合非朗伯輻射傳輸前向模型及半經(jīng)驗(yàn)核驅(qū)動模型，提出以MODIS BRDF 形狀因子先驗(yàn)數(shù)據(jù)集為驅(qū)動的地表參數(shù)確定方法，實(shí)現(xiàn)了MODIS 數(shù)據(jù)氣溶膠遙感反演。為驗(yàn)證該算法在非同源傳感器的適用性，本文基于該反演算法基本思路，以Landsat 8 OLI 為數(shù)據(jù)源，在北京地區(qū)開展AOD 反演實(shí)驗(yàn)，并對反演過程的可能誤差進(jìn)行定量分析。通過AERONET（Aerosol Robotic Network）站點(diǎn)數(shù)據(jù)及與MODIS氣溶膠產(chǎn)品交叉對比驗(yàn)證，本算法可有效提高研究區(qū)AOD 反演精度，能夠?yàn)槌鞘谐叨却髿鈴?fù)合污染防治和治理提供數(shù)據(jù)支持。

2 數(shù)據(jù)源

2.1 地面數(shù)據(jù)

AERONET作為一個(gè)全球地基氣溶膠觀測網(wǎng)絡(luò)，可實(shí)現(xiàn)高精度氣溶膠參數(shù)觀測（Holben 等，1998）。該網(wǎng)絡(luò)可提供無云條件下光譜范圍為0.34—1.06 μm的AOD 信息，同時(shí)保證數(shù)據(jù)的不確定性較低（約0.01—0.02）及高時(shí)間分辨率（約15 min）（Smirnov等，2000）。AERONET將AOD產(chǎn)品分為3級：原始Level 1.0數(shù)據(jù)、經(jīng)過云檢測Level 1.5數(shù)據(jù)及嚴(yán)格質(zhì)量控制的Level 2.0 數(shù)據(jù)。目前，最新Version 3（V3）版本在V2 基礎(chǔ)上，云檢測及數(shù)據(jù)質(zhì)量控制得到進(jìn)一步優(yōu)化（Giles 等，2019）。AERONET 已被廣泛用于AOD 衛(wèi)星反演驗(yàn)證中（Chu 等，2002；Sayer 等，2014）。本研究搜集了研究區(qū)2013 年—2018 年Beijing （BJS）、 Beijing_CAMS （BJC）、Beijing_RADI（BJR）及Xianghe（XHS）4 個(gè)站點(diǎn)V3 AOD 數(shù) 據(jù)。其 中，BJS、BJC 及XHS 站 點(diǎn) 為Level 2.0數(shù)據(jù)，BJR站點(diǎn)為Level 1.5數(shù)據(jù)。

2.2 遙感數(shù)據(jù)

Landsat 8 陸地衛(wèi)星于2013 年2 月11 日發(fā)射升空，搭載OLI（Operational Land Imager） 及TIRS（Thermal Infrared Sensor）兩個(gè)傳感器。OLI傳感器有從可見光到短波紅外9個(gè)波段，幅寬為185 km，重訪周期為16 d，空間分辨率為30 m，此外還包括一個(gè)15 m 分辨率的全色波段。研究中利用2013 年—2018年云覆蓋較小的OLI數(shù)據(jù)進(jìn)行AOD 反演實(shí)驗(yàn)，共收集到64景。MODIS傳感器搭載于Terra及Aqua衛(wèi)星上，基于該傳感器有40多種遙感參數(shù)產(chǎn)品。本文使用的MODIS 產(chǎn)品包括：AOD 產(chǎn)品（MOD04、MCD19A2）、地表反射率產(chǎn)品（MOD09A1）及地表BRDF 核驅(qū)動模型系數(shù)產(chǎn)品（MCD43A1）。其中，MOD04根據(jù)采用的反演算法不同分為DT、DB及兩者融合的DTDB AOD 產(chǎn)品，其空間分辨率均為10 km；MCD19A2 產(chǎn)品采用MAIAC 算法，空間分辨率為1 km。MOD09A1 及MCD43A1 產(chǎn)品空間分辨率均為500 m，使用該兩種產(chǎn)品構(gòu)建了地表反射率數(shù)據(jù)集及地表BRDF形狀指示因子數(shù)據(jù)集，分別實(shí)現(xiàn)HARLS 算法及本文新算法的氣溶膠反演。表1 給出了研究中使用的地面及遙感數(shù)據(jù)詳細(xì)信息。

表1 研究中使用的地面及遙感數(shù)據(jù)詳細(xì)信息Table 1 Details of ground and remote sensing data used in the study

3 原理與方法

3.1 基本原理

衛(wèi)星傳感器觀測到的大氣層頂TOA（Top Of the Atmosphere）入瞳輻亮度來自于地氣系統(tǒng)對太陽輻射的后向散射，可分為兩部分，一部分來源于大氣后向散射（程輻射項(xiàng)），另一部分則是太陽下行輻射經(jīng)地氣系統(tǒng)作用后重新散射到大氣頂層。大氣輻射傳輸模型在模擬電磁波與地氣相互作用時(shí)，通常假設(shè)地表為朗伯體。然而，真實(shí)地表為各向異性反射的非朗伯體，地表朗伯假設(shè)會影響到TOA輻射的估算，進(jìn)而影響AOD反演。為了估算地表非朗伯下TOA 反射率，Qin 等（2001）提出了考慮地表方向反射特性的輻射傳輸模型CASBIR（Correction for Anisotropic Surface BIdirectional Reflection）。該模型將輻射過程分為上行和下行分量，并且每個(gè)分量中包含直射和漫射兩部分，TOA 反射率（ρTOA）可表示為

式中，T(θs)、T(θv)為大氣下行、上行透過率矩陣，R(θs，θv，φ)為地表反射率矩陣：

由于觀測幾何及地表非朗伯特性的不同，|R|可能為正值也可能為負(fù)值。當(dāng)R的4個(gè)變量用單一反射率ρs刻畫時(shí)，|R|= 0，式（1）可簡化為朗伯假設(shè)下ρTOA計(jì)算公式，即：

式中，τ(θs)、τ(θv)為大氣下行、上行總透過率。

由式（1）可知，在計(jì)算ρTOA時(shí)需要確定與大氣有關(guān)的透過率、程輻射等大氣參數(shù)及與地表有關(guān)的反射率參數(shù)。大氣參數(shù)根據(jù)成像幾何條件，從利用輻射傳輸代碼預(yù)先構(gòu)建的反演查找表LUT（Look?Up Table）中查算獲取，LUT 構(gòu)建過程將在3.3.1 節(jié)介紹。地表參數(shù)則通過半經(jīng)驗(yàn)核驅(qū)動模型及構(gòu)建的地表BRDF 形狀因子先驗(yàn)數(shù)據(jù)集估算獲得。

3.2 地表反射率估算

3.2.1 估算模型

地表對電磁波的反射具有明顯的方向性，BRDF 可以定量地描述地表反射的各向異性特征（張虎等，2015）。根據(jù)模型構(gòu)建原理不同，BRDF可分為物理模型、經(jīng)驗(yàn)?zāi)Ｐ?、半?jīng)驗(yàn)?zāi)Ｐ图坝?jì)算機(jī)模擬模型。半經(jīng)驗(yàn)?zāi)Ｐ陀捎趨?shù)簡單，且具有一定的物理意義，已得到廣泛應(yīng)用（Lucht 等，2000），尤其是半經(jīng)驗(yàn)核驅(qū)動模型，其可表示為幾個(gè)核函數(shù)加權(quán)和的形式：

式中，λ為波長；Kvol、Kgeo分別為體散射核函數(shù)和幾何光學(xué)核函數(shù)；fiso、fvol及fgeo為核系數(shù)，分別表示各向同性散射、體散射及幾何光學(xué)散射。研究表明，羅斯表層核（Ross?Thick）和李氏稀疏互易核（LiSparse?R）核函數(shù)組合模式的擬合能力最好（Pokrovsky 和Roujean，2003），MODIS BRDF/Albedo產(chǎn)品也正是采用了該組合模式。

在已知核系數(shù)的基礎(chǔ)上，式（3）中的rdh、rhd及rhh可以通過積分獲得，Strahler 等（1999）通過大量實(shí)驗(yàn)給出了多項(xiàng)式積分?jǐn)M合式（6），表2給出了多項(xiàng)式的擬合系數(shù)。

表2 多項(xiàng)式gk = g0，k + g1，kθ2 + g2，kθ3系數(shù)取值Table 2 Coefficients for the polynomial gk = g0，k + g1，kθ2 + g2，kθ3

3.2.2 地表先驗(yàn)數(shù)據(jù)集構(gòu)建方法

由以上可知，式（4）中地表反射率估算的關(guān)鍵是核系數(shù)的確定。MCD43A1 產(chǎn)品雖然提供3 個(gè)核系數(shù)，但在反演過程中算法受到多種因素影響，包括大氣參數(shù)不確定性導(dǎo)致的反射率誤差、擬合數(shù)據(jù)的角度范圍不足及模型擬合方法無法有效抵抗噪聲等（Liang 等，2002），且受云、云陰影等影響，產(chǎn)品存在數(shù)據(jù)缺失、質(zhì)量低等問題，特別是在城市地區(qū)（Tian 等，2021b）。因此，本文利用長時(shí)間序列MCD43A1 產(chǎn)品，通過離散余弦變換的懲罰最小二乘估計(jì)DCT?PLS （Discrete Cosine Transform?Penalized Least Square regression） 時(shí) 空濾波算法進(jìn)行數(shù)據(jù)合成（Garcia 等，2010），并采用與地物類型更為相關(guān)的BRDF形狀指示因子（體散射核形狀因子：Avol=fvolfiso、幾何光學(xué)核形狀因子：Ageo=fgeofiso）進(jìn)行先驗(yàn)數(shù)據(jù)集的構(gòu)建（Tian等，2021a）。

DCT?PLS 算法的基本思想是通過計(jì)算殘差平方和及反映平滑數(shù)據(jù)粗糙度的懲罰項(xiàng)來平衡數(shù)據(jù)保真度，其對異常點(diǎn)較為敏感（Garcia 等，2010）。為降低缺失值及離群值對擬合的偏差，可使用迭代分配權(quán)重來避免，即缺失值或離群值權(quán)重低，質(zhì)量高的數(shù)據(jù)則權(quán)重高。這里參考Xiao 等（2015）在重構(gòu)地表反射率數(shù)據(jù)時(shí)采用的Bisquare加權(quán)函數(shù)進(jìn)行權(quán)重分配。

3.2.3 地表先驗(yàn)數(shù)據(jù)集構(gòu)建結(jié)果

基于上述數(shù)據(jù)濾波方法，本文使用2013 年—2017 年研究區(qū)MCD43A1 產(chǎn)品構(gòu)建了8 d 間隔、500 m分辨率地表BRDF形狀指示因子先驗(yàn)數(shù)據(jù)集。圖1 給出了第17 天及第97 天結(jié)果，可以看出，數(shù)據(jù)集在植被及城市地區(qū)均有較好的覆蓋度。由于城市結(jié)構(gòu)復(fù)雜、反射率較高，其fiso變化特征明顯高于植被和耕地地區(qū)。植被地區(qū)Avol大于城市地區(qū)，而Ageo在空間分布上整體較為均勻。利用該數(shù)據(jù)集，結(jié)合遙感影像成像幾何條件，通過式（5）和（6）即可計(jì)算獲取式（1）中各地表反射率參數(shù)，圖2 給出了一個(gè)計(jì)算結(jié)果實(shí)例。從圖中可以看出，北京及周邊地區(qū)4種反射率整體分布趨勢一致，值在0.04—0.10 之間，各反射率之間存在明顯差異，其中rhh大于其他反射率。

圖1 北京及周邊地區(qū)BRDF形狀指示因子先驗(yàn)數(shù)據(jù)集Fig.1 A priori dataset of BRDF shape factor parameters over Beijing and surrounding areas

圖2 北京及周邊地區(qū)4種地表反射率參數(shù)估算結(jié)果Fig.2 Estimation results of four surface reflectances over Beijing and surrounding areas

3.2.4 先驗(yàn)數(shù)據(jù)集擬合能力評估

為評估本文構(gòu)建的先驗(yàn)數(shù)據(jù)集擬合地表反射率的能力，我們對研究區(qū)3 種典型地物類型（耕地、人工地表及自然植被）進(jìn)行了反射率擬合實(shí)驗(yàn)，結(jié)果如圖3。其中，估算地表反射率采用構(gòu)建的先驗(yàn)數(shù)據(jù)集，結(jié)合半經(jīng)驗(yàn)核驅(qū)動模型計(jì)算獲得，觀測地表反射率數(shù)據(jù)則為MODIS 逐日地表反射率產(chǎn)品（MOD09GA）。從圖3中可以看出，估算的地表反射率與觀測值具有較好的一致性，尤其在耕地及自然植被地表，相關(guān)系數(shù)R大于0.985，RMSE值小于0.005；在人工地表，反射率估算結(jié)果平均絕對誤差MAE為0.005。

圖3 BRDF形狀指示因子先驗(yàn)數(shù)據(jù)集擬合地表反射率與MODIS反射率產(chǎn)品對比散點(diǎn)圖Fig.3 Scatter plot of surface reflectance fitted from a priori dataset of BRDF shape factor parameters compared to MODIS reflectance product

3.3 氣溶膠光學(xué)厚度反演

3.3.1 查找表構(gòu)建

氣溶膠及大氣分子共同決定了式（1）中與大氣有關(guān)的參數(shù)，大氣氣溶膠的組成隨時(shí)間和地點(diǎn)變化而不同，氣溶膠遙感另一個(gè)關(guān)鍵是氣溶膠類型的確定（李正強(qiáng)等，2019）。MODIS氣溶膠反演算法根據(jù)季節(jié)和位置從預(yù)先建立的全球氣溶膠類型數(shù)據(jù)庫中選擇合適的氣溶膠類型（Omar 等，2005）。研究發(fā)現(xiàn)，北京及周邊地區(qū)氣溶膠狀況可以用大陸型氣溶膠描述（Sun 等，2016；田信鵬等，2018）。本文利用MODTRAN 5 輻射傳輸代碼構(gòu)建通用大氣參數(shù)數(shù)據(jù)集，該數(shù)據(jù)集包含不同觀測 幾 何 及AOD 條 件 下0.3—2.5 μm 光 譜 范 圍、15 cm?1光譜分辨率的6 個(gè)大氣輸出參數(shù)，如表3 所示?；诖藬?shù)據(jù)集，結(jié)合傳感器光譜響應(yīng)函數(shù)即可快速構(gòu)建某特定傳感器AOD反演查找表。

表3 基于MODTRAN構(gòu)建的通用大氣參數(shù)數(shù)據(jù)集Table 3 General atmospheric parameters dataset based on MODTRAN

3.3.2 像元篩選

氣溶膠遙感反演的成功不僅與地表反射率的精確確定有關(guān)，在很大程度上也依賴于不合適像元的有效剔除，例如云、云陰影、冰雪及內(nèi)陸水體等（Levy 等，2013）。Zhu 等（2015） 針對Landsat 8 OLI 傳 感 器 提 出 的Fmask （Function of mask）算法可實(shí)現(xiàn)云、云陰影、雪及內(nèi)陸水體檢測。本文在反演AOD 前，首先基于該算法進(jìn)行可反演區(qū)域的像元篩選。

3.3.3 光譜轉(zhuǎn)換

不同傳感器間光譜響應(yīng)條件存在差異，導(dǎo)致具有相同反射特性的地物在不同傳感器表現(xiàn)為不同特征。在利用由先驗(yàn)數(shù)據(jù)集支持非同源其它傳感器AOD 反演時(shí)，需將先驗(yàn)數(shù)據(jù)集的反射率轉(zhuǎn)化為待反演傳感器對應(yīng)的反射率，消除兩傳感器之間的光譜差異（Tian 等，2018）。為將MODIS 反射率轉(zhuǎn)化為OLI 波段反射率，本文選用《定量遙感：理念與算法》（梁順林等，2013）所附光盤中波譜數(shù)據(jù)及北京師范大學(xué)“中國典型地物標(biāo)準(zhǔn)波普數(shù)據(jù)庫”植被和土壤波譜數(shù)據(jù)作為樣本，擬合獲取OLI 與MODIS 藍(lán)關(guān)波段反射率轉(zhuǎn)換關(guān)系，結(jié)果為：

ρOLIs= 1.047 ×ρMODISs+ 0.001。

圖4給出了本文AOD反演算法的具體流程。

圖4 反演算法流程圖Fig.4 Flowchart of AOD retrieval in this study

4 結(jié)果與討論

4.1 空間分布對比

圖5 對比展示了不同算法獲取的研究區(qū)AOD空間分布情況，包括本文算法（AOD_BRDF）、地表朗伯假設(shè)（WSA 反射率支持） 反演（AOD_WSA）、HARLS 算 法（AOD_HARLS） 及MODIS DT、DB、DTDB、MAIAC 產(chǎn)品。整體來看，各算法獲取的AOD 在空間分布上具有較好的一致性，均表現(xiàn)為城市地區(qū)高、植被地區(qū)低的特點(diǎn)。其中，AOD_BRDF 整體小于AOD_WSA 及AOD_HARLS，MODIS DT及DB產(chǎn)品較其他反演結(jié)果偏大。

圖5 不同算法反演獲取的AOD空間分布對比圖（2013年5月12日）Fig.5 Comparisons of spatial distributions of AOD obtained by different algorithms（May 12th，2013）

4.2 精度驗(yàn)證

地面測量可提供站點(diǎn)周圍高質(zhì)量、高頻率氣溶膠數(shù)據(jù)，而衛(wèi)星遙感則提供過境時(shí)刻、一定空間范圍氣溶膠綜合狀況。利用地面觀測資料驗(yàn)證遙感產(chǎn)品時(shí)，需對二者進(jìn)行時(shí)空匹配（Yang 等，2014）。一般采取的空間匹配策略為選取站點(diǎn)周圍一定大小窗口內(nèi)遙感產(chǎn)品平均值作為最終反演值，時(shí)間匹配策略為將衛(wèi)星過境前后一定時(shí)間段的地面測量平均值作為最終測量值（Levy 等，2010）。因此，在后續(xù)驗(yàn)證過程中，本文對反演結(jié)果進(jìn)行5×5 空間平均處理，對AERONET 站點(diǎn)數(shù)據(jù)進(jìn)行衛(wèi)星過境前后30 min 時(shí)間平均處理。為定量描述算法反演精度，采用的評價(jià)指標(biāo)包括：匹配點(diǎn)個(gè)數(shù)（N），相關(guān)系數(shù)（R），均方根誤差（RMSE）及定義 為± （0.05+0.15×AODAERONET） 的 期 望 誤 差EE（Expected Error）。此外，為分析算法的高估和低估情況，將在EE 誤差線以內(nèi)（Within）的點(diǎn)對百分比稱為PWE，在EE 誤差線以上（Above）和以下（Below）的點(diǎn)對占比分別稱為PAE和PBE。

4.2.1 AERONET驗(yàn)證

由于BJS、BJC及BJR站點(diǎn)地表類型均為城市，且3 個(gè)站點(diǎn)距離較近。受云污染及數(shù)據(jù)量限制，BJS、BJC 及BJR 站點(diǎn)N分別為39、41、37，XHS站點(diǎn)N= 43。鑒于此，本文將BJS、BJC 及BJR 站點(diǎn)作為城市站點(diǎn)對其進(jìn)行整體評價(jià)，XHS 站點(diǎn)作為植被站點(diǎn)進(jìn)行評價(jià)分析。

圖6展示了城市站點(diǎn)不同反演算法反演值與地面實(shí)測值對比結(jié)果，其中紅色實(shí)線為擬合線，黑色實(shí)線為1∶1線，黑色虛線為EE 誤差線；從中可以看出，AOD_HARLS（圖6（c））各精度指標(biāo)（PWE = 73.5%，R= 0.978，RMSE = 0.164） 比AOD_WSA（圖6（b））稍有提高（PWE = 65.8%，R= 0.977，RMSE = 0.173），兩者PAE 值分別為22.2%、27.4%，擬合線斜率分別為1.19、1.21，表明對AOD 均存在高估。圖6（a）為新算法反演值與實(shí)測值對比結(jié)果，可以看出，兩者具有較高的一致性（R= 0.986，RMSE = 0.102），PWE 值達(dá)到84.6%，較AOD_WSA、AOD_HARLS 分別提高了18.8%、11.1%。AOD_BRDF 結(jié)果改善主要體現(xiàn)在對高估的有效抑制，其PAE =14.5%，比AOD_WSA、AOD_HARLS分別低了12.9%、7.7%。

圖6 城市站點(diǎn)（BJS、BJC及BJR）反演結(jié)果與AERONET對比散點(diǎn)圖Fig.6 Scatter plots of AOD retrievals against AERONET measurements from urban stations of BJS，BJC and BJR

圖7 為XHS 植被站點(diǎn)不同算法獲取的AOD 反演值與地面實(shí)測值對比散點(diǎn)圖，從中可以看出，AOD_BRDF 精度最高（PWE = 86.0%，R= 0.991，RMSE=0.086），PWE 較AOD_WSA 及AOD_HARLS分別提高了20.9%和9.3%；其次為AOD_HARLS（PWE=76.7%，R=0.985，RMSE=0.119），AOD_WSA反演精度最低（PWE=65.1%，R=0.976，RMSE=0.165）。與城市站點(diǎn)一樣，AOD_BRDF精度的提高也體現(xiàn)在對AOD 高估的抑制，其PAE 值比AOD_WSA、AOD_HARLS分別低了18.6%、9.3%。

圖7 XHS植被站點(diǎn)反演結(jié)果與AERONET對比散點(diǎn)圖Fig.7 Scatter plots of AOD retrievals against AERONET measurements from vegetation station of XHS

以上分析表明，在城市及植被地區(qū)，本文提出的反演算法均可有效改善地表朗伯假設(shè)（即使用WSA 反射率支持）及HARLS 算法的高估，提高AOD遙感反演精度。

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

圖8（a）—8（d）給出了AOD_BRDF 與MODIS DT、DB、DTDB 及MAIAC 產(chǎn)品在城市站點(diǎn)的精度對比，其中紅色點(diǎn)為新算法反演結(jié)果，藍(lán)色點(diǎn)為相應(yīng)的MODIS 產(chǎn)品。從圖8 中可以看出，DT、DB、DTDB 及MAIAC 算法N分別為47、79、61 及92，表明各算法中，MAIAC 算法在城市地區(qū)適用性最好，其次為DB 算法，DT 算法適用性最差，但均要小于本文AOD_BRDF 結(jié)果（N= 117，見圖6）。與DT算法相比，AOD_BRDF 與實(shí)測值具有更好的一致性：R由0.975 提高到0.996，RMSE 值由0.154 降為0.117。新算法有效降低了DT 算法的高估，PAE 由57.4%減小到10.6%，PWE 提高了46.8%；DB 與DTDB 算法 的PWE 分別為70.9%、62.3%，與這兩種產(chǎn)品相比，AOD_BRDF 同樣有不同幅度精度提升，PWE 分別提高了13.9% 及14.7%。新算法精度的提高同樣表現(xiàn)在對高估現(xiàn)象的改善，PAE 分別降低了5.1%及14.7%。MODIS各算法中，MAIAC 算法在城市站點(diǎn)具有最高的精度（PWE=79.3%），本文新算法與之相比PWE 提高了4.4%，精度也有小幅度提升。

圖8 城市站點(diǎn)本文算法與MODIS氣溶膠產(chǎn)品對比圖Fig.8 Scatter plots of the new method vs.MODIS AOD retrievals against AERONET measurements from urban stations of BJS，BJC and BJR

圖9 展示了本文算法與MODIS AOD 產(chǎn)品在XHS植被站點(diǎn)的交叉對比。通過與城市站點(diǎn)對比可以看出，DT算法在植被地區(qū)有了較大改善（PWE=68.4%）（圖9（a））。DT、DB 及DTDB 算法在XHS站點(diǎn)對AOD同樣存在高估現(xiàn)象，PAE分別為31.6%、41.9%及48.3%。與DT、DB及DTDB相比，本文新算法在該站點(diǎn)精度有不同幅度提升：PAE分別降低了21.1%、22.5%及27.6%，PWE 提高超過21.1%。新算法與MAIAC相比，PWE均為82.9%，但具有更高的R值（0.992）及更小的RMSE（0.085），表明新算法與實(shí)測值有更好的一致性。

圖9 XHS植被站點(diǎn)本文算法與MODIS氣溶膠產(chǎn)品對比圖Fig.9 Scatter plots of the new method vs.MODIS AOD retrievals against AERONET measurements from vegetation station of XHS

通過以上分析表明，本文反演算法比WSA 反射率支持及HARLS算法反演結(jié)果具有更高的精度。由于AOD_BRDF、AOD_WSA 及AOD_HARLS 反演中使用的LUT 相同，因此，新算法反演精度的提高是地表參數(shù)估算優(yōu)化的結(jié)果。此外，新算法在城市及植被站點(diǎn)均表現(xiàn)出比MODIS 產(chǎn)品更好的適用性，可以有效地改善當(dāng)前氣溶膠遙感反演在研究區(qū)的高估。

4.3 不確定性分析

4.3.1 地表反射率影響

本文新算法、WSA 反射率支持及HARLS 算法采用不同地表反射率估算策略進(jìn)行AOD 反演，圖10 給出了BJS 及XHS 站點(diǎn)沿主平面方向，觀測天頂角θv= 0°時(shí)的方向?方向反射率（rdd）、半球?半球反射率（rhh，WSA）及MOD09A1 產(chǎn)品合成反射率（rmod09）隨太陽天頂角（θs）的變化情況。rhh是通過構(gòu)建的地表BRDF形狀因子先驗(yàn)數(shù)據(jù)集利用式（6）計(jì)算獲得，rmod09為HARLS算法使用的地表反射率。從圖中可以看出，BJS、XHS 站點(diǎn)的rdd在主平面方向均隨θs遠(yuǎn)離“熱點(diǎn)”角度（θs=θv）減小，在“熱點(diǎn)”方向達(dá)到最大值0.064、0.048。隨著θs的變化，兩站點(diǎn)rdd值變化幅度超過0.030。Landsat 8 影像行/列號為163/47 時(shí)，OLI 傳感器成像θs的范圍約為20°—65°（圖10 中虛線）。在該角度內(nèi)，兩站點(diǎn)rdd值變化范圍為0.025。因此，在AOD 反演中若不考慮地表非朗伯特性對地表反射率估算的影響勢必會帶來較大的誤差。

圖10 沿主平面方向不同地表反射率隨太陽天頂角的變化（θv=0°）Fig.10 Variation of surface reflectances with solar zenith angle in the principal plane for θv of 0°

WSA反射率支持的AOD反演及HARLS算法在反演過程中使用的rhh和rmod09均為固定值，且BJS及XHS 兩站點(diǎn)的rhh值小于rmod09。因此，圖6 城市站點(diǎn)和圖7 植被站點(diǎn)中AOD_WSA 比AOD_HARLS有更嚴(yán)重高估主要原因來自于地表反射率因素。以上分析表明，由于城市和植被下墊面方向反射特性，導(dǎo)致不同成像幾何條件下地表反射率具有較大差異。在AOD 反演過程中如果使用單一角度或合成反射率來刻畫地表信息，將會給反演帶來較大誤差。

4.3.2 空間尺度轉(zhuǎn)換影響

由于地表異質(zhì)的普遍性，單個(gè)遙感像元往往包含多種地表覆蓋類型（陳晉等，2016）。本文利用Landsat 8 OLI 30 m 空間分辨率數(shù)據(jù)反演氣溶膠，算法中地氣分離策略是基于500 m 分辨率先驗(yàn)數(shù)據(jù)集。由于OLI數(shù)據(jù)與先驗(yàn)數(shù)據(jù)集存在空間尺度不一致問題，反演過程中將OLI高空間分辨率數(shù)據(jù)升尺度到與先驗(yàn)數(shù)據(jù)集相同的500 m 分辨率，重采樣會使遙感數(shù)據(jù)的表觀反射率發(fā)生變化，進(jìn)而影響AOD反演。

為定量分析空間升尺度所帶來的反演不確定性，我們分析了研究區(qū)3種典型地物（耕地、人工地表和自然植被）空間尺度轉(zhuǎn)化對TOA 反射率重采樣的影響，如圖11（a）所示。從中可以看出，3 種地物由30 m 升尺度到500 m 分辨率，其TOA反射率絕對差值分別為0.0029、0.0017 和0.0013。圖11 （b） 為基于輻射傳輸模型模擬獲取的不同AOD 及地表反射率條件下的表觀反射率。當(dāng)AOD<1.5 時(shí)，耕地、人工地表和自然植被擬合斜率分別為0.040、0.034 和0.052。因此，上述表觀反射率絕對差值對應(yīng)的AOD 誤差約為0.073、0.050和0.025。可見耕地和人工地表在升尺度過程導(dǎo)致的AOD 誤差最大，這可能是由于北京地區(qū)該地表類型空間異質(zhì)性較自然植被更高，升尺度所帶來的表觀反射率差異更大造成的。

圖11 空間尺度轉(zhuǎn)換對AOD反演精度影響模擬分析Fig.11 Simulation analysis of spatial scale conversion effect on AOD inversion accuracy

5 結(jié) 論

本文結(jié)合非朗伯大氣輻射傳輸前向模型及半經(jīng)驗(yàn)核驅(qū)動模型，提出以地表先驗(yàn)數(shù)據(jù)集為驅(qū)動的地表參數(shù)確定方法，實(shí)現(xiàn)了考慮地表BRDF效應(yīng)的氣溶膠遙感反演。利用Landsat 8 OLI 數(shù)據(jù)在北京地區(qū)進(jìn)行了反演實(shí)驗(yàn)，得到主要結(jié)論如下：

（1）利用長時(shí)間序列MODIS BRDF/Albedo 產(chǎn)品，基于離散余弦變換的懲罰最小二乘估計(jì)（DCT?PLS）時(shí)空濾波算法構(gòu)建了8 d 間隔、500 m分辨率地表BRDF形狀指示因子先驗(yàn)數(shù)據(jù)集。

（2）利用MODTRAN 5 輻射傳輸代碼構(gòu)建了通用大氣參數(shù)數(shù)據(jù)集，基于該數(shù)據(jù)集，結(jié)合傳感器光譜響應(yīng)函數(shù)即可快速構(gòu)建查找表。

（3）通過分析非朗伯輻射傳輸前向模型及半經(jīng)驗(yàn)核驅(qū)動模型特點(diǎn)，實(shí)現(xiàn)了兩模型的一體化耦合，提出了考慮地表BRDF 效應(yīng)的AOD 反演算法，突破傳統(tǒng)算法的地表朗伯假設(shè)。

（4）基于本文算法實(shí)現(xiàn)了Landsat 8 OLI 數(shù)據(jù)在北京地區(qū)AOD 反演實(shí)驗(yàn)。通過與朗伯假設(shè)反演及HARLS 算法反演結(jié)果的交叉對比，結(jié)果表明，新算法在城市/植被地區(qū)誤差線內(nèi)比重達(dá)到了84.6%/86.0%，分別較前兩種算法提高了18.8%/20.9%和11.1%/9.3%。由于3 種算法反演中采用相同的LUT，因此，新算法精度提高是地表參數(shù)優(yōu)化的結(jié)果。

（5）通過與MODIS DT、DB、DTDB 及MAIAC氣溶膠產(chǎn)品交叉對比，新算法在城市及植被站點(diǎn)均表現(xiàn)出與實(shí)測值更好的一致性及反演精度。

（6）誤差分析表明，由于地表方向反射會導(dǎo)致不同成像幾何條件下地表反射率差異較大，在AOD 反演中若忽略地表非朗伯特性，使用單一角度或合成反射率來刻畫地表信息，會帶來較大不確定性。