中國陸表穩(wěn)定目標TOA反射率模型構(gòu)建方法及驗證
——以格爾木沙地場景為例

宋培蘭，馬靈玲，趙永光，王寧，李婉，韓啟金，劉耀開，姚微源，張貝貝，任璐，牛沂芳

1.中國科學(xué)院空天信息創(chuàng)新研究院 中國科學(xué)院定量遙感信息技術(shù)重點實驗室,北京 100094;

2.中國科學(xué)院大學(xué) 光電學(xué)院,北京 100049;

3.中國資源衛(wèi)星應(yīng)用中心,北京 100094

1 引言

偽不變定標場PICS（Pseudo?Invariant Calibration Sites）定標通過嚴格的條件篩選地表/大氣特性極其穩(wěn)定的地球目標區(qū)域，并構(gòu)建高精度PICS大氣層頂TOA（Top of Atmosphere）反射率模型作為其輻射參考，實現(xiàn)載荷輻射性能長時序監(jiān)測和跟蹤。自1984 年定標與真實性檢驗工作組WGCV（Working Group on Calibration &Validation）成立以來，國際衛(wèi)星對地觀測委員會CEOS（Committee of Earth Observation Satellite）就一直致力于全球范圍內(nèi)輻射定標場的優(yōu)選，并于2008 年確定了Algeria 3、Algeria 5、Libya 1、Libya 4、Mauritania 1 及Mauritania 2共6個PICS（Lachérade 等，2013），已被廣泛用于Terra/MODIS、Landsat 7 ETM+、FY?3A MERSI、EO?1/Hyperion、OCO?2 等系列傳感器長期穩(wěn)定性評估（Chander 等，2010；Kim 等，2014；Neigh 等，2016；Yu 等，2020）。然而，PICS 有著嚴格的約束條件，例如場地面積大于20 km×20 km、云覆蓋率<20%（1 年）、空間均勻性<3%、時間穩(wěn)定性<4%等。截至目前，Bacour等（2019）最新甄選的PICS均分布于國外北非撒哈拉沙漠地區(qū)，處于北回線附近，屬于典型的熱帶沙漠氣候，地表裸露、終年干旱炎熱、四季變化不明顯。在針對國外PICS 場地的TOA 反射率模型構(gòu)建研究中，Helder 等（2013）使用Terra/MODIS 和EO?1/Hyperion 建立了一種Libya4 沙漠目標簡化的絕對輻射定標模型，該模型只考慮太陽天頂角變化對TOA 反射率的影響，可計算近似星下點觀測條件（±7.5°）時目標上空的TOA 反射率。Mishra等（2014）進一步引入觀測天頂角的影響，可實現(xiàn)非星下點觀測、不同大氣條件下目標上空TOA反射率的模擬計算，與Aqua/MODIS、Landsat 8 OLI、UK?2/DMC、ENVISAT/MERSI 等衛(wèi)星觀測數(shù)據(jù)的相對誤差標準差小于3%。Raut 等（2019）將該模型推廣到Egypt 1、Libya 1、Niger 1、Niger 2等場地，建立的通用TOA 反射率模型在Libya 1 的不確定度為1.07%—4.02%。由于PICS穩(wěn)定的地表/大氣特性，使得僅使用太陽天頂角和觀測天頂角的組合就能夠表征PICS 場地TOA 反射率的規(guī)律性變化，并在應(yīng)用中達到較好效果。

由于中國陸地衛(wèi)星成像資源多在國內(nèi)，因此國內(nèi)學(xué)者更多關(guān)注中國境內(nèi)的陸表穩(wěn)定目標甄選，何靈莉（2020）以時間穩(wěn)定性<5%（16 年），空間均勻性<5%為選場約束，選場結(jié)果主要集中在內(nèi)蒙古以及西北區(qū)域。受溫帶大陸性氣候及高原山地氣候影響，這些場地的地表、大氣特性不如PICS 場地穩(wěn)定。因此，建立于地表/大氣特性穩(wěn)定性極高基礎(chǔ)之上的PICS 場地TOA 反射率模型難以準確表征國內(nèi)陸表穩(wěn)定目標由角度、大氣以及季節(jié)等因素引起的TOA 反射率規(guī)律性變化。為滿足中國多系列陸地衛(wèi)星的定標需求，有必要針對國內(nèi)目標特性相對穩(wěn)定的場地開展場地特性研究，發(fā)展適用于國內(nèi)穩(wěn)定目標場地的TOA 反射率模型，開展高頻次的衛(wèi)星在軌輻射定標和性能監(jiān)測，將是對國產(chǎn)衛(wèi)星業(yè)務(wù)化場地定標的有益補充，對提高其輻射定標頻次和穩(wěn)定性具有重要意義。近幾年，一些學(xué)者在國內(nèi)穩(wěn)定目標TOA 反射率建模研究中，借鑒地表BRDF 模型（趙春艷 等，2019）建立的經(jīng)驗，基于中國敦煌輻射定標場和格爾木場地一年中晴天MODIS TOA 反射率和角度信息建立TOA 反射率半經(jīng)驗核驅(qū)動模型（Liu 等，2019），但使用兩個場地分別得出的定標系數(shù)一致性較差，在綠波段（B4）的相對差異高達7.70%。Wang 等（2017）考慮下墊面地物類型選擇雪面BRDF 模型對Kunlun1的TOA反射率進行建模，未考慮大氣季節(jié)性變化對TOA 反射率的影響，導(dǎo)致模型在夏季出現(xiàn)較大偏差。此外，王詩圣等（2020）基于6S模型分別建立了TOA 反射率隨觀測天頂角和相對方位角變化的經(jīng)驗公式，但未考慮太陽角度影響?？偟膩碚f，國內(nèi)穩(wěn)定目標TOA 反射率模型構(gòu)建多是從描述地表二向反射特性的半經(jīng)驗?zāi)Ｐ统霭l(fā)，均忽略了國內(nèi)場地大氣季節(jié)性變化對TOA 反射率的影響，難以表征國內(nèi)穩(wěn)定目標由角度、大氣以及季節(jié)等多因素綜合影響引起的TOA 反射率規(guī)律性變化。

本文提出一種考慮“太陽—目標—傳感器”觀測幾何及地表/大氣季節(jié)性變化的TOA 反射率模型構(gòu)建方法，在放寬場地目標空間均勻性、地表/大氣穩(wěn)定性以及BRDF特性等約束條件下，選取了中國青海省柴達木盆地周邊5 km×5 km 的格爾木沙地穩(wěn)定目標，以MODIS 和ECMWF 數(shù)據(jù)產(chǎn)品作為穩(wěn)定目標的TOA 反射率、大氣狀況先驗知識庫，分析觀測幾何以及大氣參數(shù)對TOA 反射率的影響規(guī)律，構(gòu)建陸表穩(wěn)定目標TOA 反射率模型，并利用高精度的Landsat 8 OLI、Sentinel?2A/B MSI 觀測數(shù)據(jù)進行了定標結(jié)果分析，以驗證所發(fā)展模型的有效性。

2 研究區(qū)域及數(shù)據(jù)

2.1 研究區(qū)域概況

借鑒國外PICS 目標甄選標準，并參考國內(nèi)穩(wěn)定目標選取的經(jīng)驗和規(guī)則，本文首先制定了國內(nèi)陸表穩(wěn)定目標的約束條件。其中，場地面積從PICS 的20 km×20 km 放寬至3 km×3 km；年度云覆蓋率小于60%以便獲取足夠的晴空影像，晴空條件下地表反射率應(yīng)在VIS>0.15，NIR>0.3；當(dāng)太陽及傳感器天頂角小于65°時，角度效應(yīng)小于60%；空間均勻性應(yīng)優(yōu)于3%，時間穩(wěn)定性<10%（10 年以上）。

本文基于2020 年3—10 月Sentinel?2A2B 地表反射率圖像，以3 km × 3 km 的滑動窗口評價了過境中國青海省格爾木市附近區(qū)域的19 景晴空影像的空間均勻性，以3%為閾值篩選出5 km×5 km 的格爾木穩(wěn)定目標。其下墊面類型為沙地，如圖1所示。圖2為在研究區(qū)實地拍攝的沙地圖像，其顆粒較粗。2019年8月9日及14日分別實測了格爾木沙地穩(wěn)定目標的反射率光譜，如圖3所示，約20%—40%。經(jīng)2010年—2020年地表反射率產(chǎn)品MCD43A4分析，研究區(qū)域在可見—近紅外波段（B3、B4、B1、B2）的時間穩(wěn)定性分別為8.37%、6.70%、5.26%、5.07%。距格爾木市約50 km，地理坐標為：36.41°N—36.46°N，94.20°E—94.25°E，總面積約5 km×5 km，中心海拔高度約2.84 km。海拔高、大氣干潔穩(wěn)定，年平均晴空MODIS 影像數(shù)量約為150天。

圖1 格爾木沙地的Sentinel?2B L2A影像（2020?05?01）Fig.1 The Sentinel?2B L2A image of Golmud Desert

圖2 實地拍攝的格爾木沙地Fig.2 The image of Golmud Desert taken at the scene

圖3 格爾木沙地穩(wěn)定目標反射率光譜Fig.3 Spectral reflectance of stable target in Golmud Desert

2.2 數(shù)據(jù)介紹

穩(wěn)定目標TOA 反射率模型構(gòu)建依賴于包含衛(wèi)星TOA 反射率、大氣參數(shù)等在內(nèi)的場地特性數(shù)據(jù)，因此，針對本文研究區(qū)域，基于衛(wèi)星數(shù)據(jù)、產(chǎn)品以及再分析資料搜集了建模所需的陸表穩(wěn)定目標特性數(shù)據(jù)集，主要包括：

（1）Aqua/MODIS、Sentinel?2A/B MSI 以 及Landsat 8 OLI 衛(wèi)星數(shù)據(jù)。Aqua/MODIS 具有多角度觀測能力，利用其寬視場角帶來的軌道重疊獲取多時相數(shù)據(jù)，可累積形成多角度數(shù)據(jù)集（閻廣建 等，2021）。本研究選擇2010 年—2020 年Aqua衛(wèi)星過境格爾木沙地目標的Collection 6.0 版本MYD021KM 輻射定標產(chǎn)品、MYD03 地理定位產(chǎn)品作為研究數(shù)據(jù)，篩選出無云無雪晴空影像1232景，提取其太陽反射波段TOA 反射率及觀測幾何參數(shù)構(gòu)建模型構(gòu)建所需的穩(wěn)定目標數(shù)據(jù)集。選擇2018 年—2021 年30 景Sentinel?2A/B MSI影像及2013 年—2020年27景Landsat 8 OLI影像作為驗證數(shù)據(jù)。圖4 展示了這3 個衛(wèi)星載荷對應(yīng)波段的光譜響應(yīng)曲線，其中黃色線表示的是Sentinel?2B 衛(wèi)星上搭載的MSI的光譜響應(yīng)曲線，Sentinel?2A/MSI光譜響應(yīng)曲線與其相似。

圖4 Aqua/MODIS、Sentinel?2B/MSI及Landsat 8 OLI相匹配的波段Fig.4 The matching bands of Aqua/MODIS，Sentinel?2B/MSI and Landsat 8 OLI

（2）再分析資料。本文使用歐洲中期天氣預(yù)報中心ECMWF（European Center for Medium?Range Weather Forecasts）再分析資料產(chǎn)品，從空間分辨率為0.25°× 0.25°的ERA5 數(shù)據(jù)庫以及空間分辨率為0.75°× 0.75°的EAC4 數(shù)據(jù)庫中提取2010 年—2020年Aqua/MODIS 衛(wèi)星過境格爾木沙地目標時的陸表降水量、降水類型、水汽含量、氣溶膠光學(xué)厚度、風(fēng)速、臭氧等大氣參數(shù)。

3 方法及原理

3.1 基本思路

本文提出一種適用于國內(nèi)陸表穩(wěn)定目標的TOA反射率模型建模方法，該模型引入觀測幾何參數(shù)及大氣參數(shù)，能夠較好的表征觀測幾何、大氣條件引起的場地目標TOA 反射率變化，基本思路如圖5所示。

圖5 TOA反射率模型構(gòu)建方法Fig.5 The method of constructing TOA reflectance model

（1）考慮“太陽—目標—傳感器”觀測幾何、地表/大氣特性參數(shù)等因素，分析其對場地目標TOA 反射率的影響。包括：分析TOA 反射率隨觀測角度或太陽角度的變化規(guī)律，對TOA 反射率隨太陽及傳感器觀測角度變化進行建模；同時分析降水量、水汽含量、臭氧含量、氣溶膠光學(xué)厚度等因素的影響，對大氣季節(jié)性影響進行建模。

（2）制定包含觀測幾何、大氣參數(shù)的模型適用性約束條件，并構(gòu)建陸表穩(wěn)定目標TOA 反射率模型。主要考慮太陽、觀測角度對大氣路徑影響，并考慮風(fēng)沙、降雨、降雪等不滿足輻射定標條件的天氣。當(dāng)滿足觀測幾何、大氣參數(shù)約束條件時，以陸表穩(wěn)定目標大氣及場地特性為輻射參考，耦合觀測幾何及大氣參數(shù)，建立的TOA反射率模型為

式中，ρTOA是TOA 反射率，θs、φs、θv、φv、WVC、O3、及AOD是太陽觀測幾何、傳感器觀測幾何、水汽含量、臭氧含量、氣溶膠光學(xué)厚度等影響因素。

3.2 觀測幾何影響建模

3.2.1 觀測天頂角對TOA反射率的影響

為分析觀測天頂角引起的TOA 反射率規(guī)律性變化，本文利用了2010 年—2020 年Aqua/MODIS的TOA 反射率時間序列數(shù)據(jù)。其中，Aqua/MODIS過境格爾木沙地目標時觀測天頂角覆蓋范圍為：0°—65°，受衛(wèi)星軌道周期性變化影響，衛(wèi)星觀測角度也呈現(xiàn)規(guī)律性趨勢，主要集中在不連續(xù)的17 組數(shù)據(jù)中。圖6 展示了格爾木沙地穩(wěn)定目標TOA 反射率隨觀測天頂角的變化規(guī)律，其中衛(wèi)星和太陽同向時觀測天頂角記為正，對向觀測時記為負。從圖中展示的數(shù)據(jù)可以看出，隨著衛(wèi)星逐漸偏離垂直觀測，反射率有逐漸增大的趨勢。在不考慮太陽觀測幾何的情況下，為建立TOA 反射率與觀測天頂角θv的關(guān)系模型，使用余弦函數(shù)對不同波段的TOA 反射率數(shù)據(jù)進行擬合，擬合結(jié)果如圖6 所示。其第一項常數(shù)項不為0，表征衛(wèi)星和太陽同向觀測與對向觀測時TOA 反射率的不對稱現(xiàn)象；第二項常數(shù)項反映了對應(yīng)波段TOA 反射率均值大??；一次項表示觀測天頂角θv對各波段TOA 反射率的影響高低。對于Band3波段，θv的變化對TOA 反射率的影響最為顯著，對應(yīng)一次項絕對值最大。對于Band7波段，相應(yīng)的一次項絕對值數(shù)相對較小，但是θv對TOA反射率仍有影響。

圖6 AQUA/MODIS傳感器過境格爾木沙地時TOA反射率隨觀測天頂角的變化規(guī)律Fig.6 The TOA reflectance changes with the view zenith angle when AQUA/MODIS sensor passes through the Golmud Desert

隨著θv的變化，輻射能量穿透大氣的光學(xué)路徑長度也隨之變化。以AQUA 衛(wèi)星為例，分析了不同θv下衛(wèi)星星下點與目標場地的距離，結(jié)果如圖7 所示。當(dāng)θv=35°度時，星下點與觀測目標之間的距離約493.65 km，因此，本文通過約束θv以限制大氣路徑長度。

圖7 MODIS傳感器下視點距離觀測目標的距離Fig.7 The distance from lower viewpoint to observation target

3.2.2 太陽天頂角對TOA反射率的影響

考慮太陽天頂角變化對穩(wěn)定目標TOA 反射率的影響，本文采用控制變量法，分析并建立TOA反射率與太陽天頂角θs的參數(shù)化關(guān)系，也就是說，在限制θv及傳感器觀測方位角φv變化范圍的基礎(chǔ)上，分析TOA反射率隨θs的變化規(guī)律。

本文使用的MODIS 數(shù)據(jù)過境時間為：UTC6：00—7：40，θs覆蓋范圍為0°—65°。依據(jù)θv的數(shù)值變化，將這些數(shù)據(jù)劃分為17 組。同一組數(shù)據(jù)中θv的增量不超過2°。不同組數(shù)據(jù)間φv變化的增量為0.23°—2.60°，僅當(dāng)?3°≤θv≤?1°時，φv的變化范圍為77.10°—97.60°，增量為20.50°。因此，每組數(shù)據(jù)對應(yīng)的傳感器觀測姿態(tài)變化較小，可以假定TOA 反射率主要受太陽天頂角影響。圖8展示了其中6 組特定觀測天頂角下TOA 反射率（MODIS B7波段）隨θs的變化規(guī)律。綜合分析，當(dāng)?64°≤θv≤?41°、47°≤θv≤64°時，TOA 反射率隨θs的增大而明顯增大；當(dāng)?35°≤θv≤?24°、40°≤θv≤42°時，TOA 反射率隨θs而增加的趨勢減緩；當(dāng)?15°≤θv≤33°時，TOA 反射率隨θs的增大而明顯減小。結(jié)合3.2.1 節(jié)的分析可知，θv=35°近似同向觀測場景下數(shù)據(jù)趨勢變化的臨界值。

圖8 AQUA/MODIS傳感器過境格爾木沙地時B7的TOA反射率隨太陽天頂角的變化Fig.8 The TOA reflectance in B7 varies with the solar zenith angle when AQUA/MODIS sensor passes through the Golmud Desert

為建立TOA 反射率與θs的關(guān)系模型，使用θs余弦的一次函數(shù)對不同波段的TOA 反射率數(shù)據(jù)進行擬合。圖8 給出了B7 波段的擬合結(jié)果，在限定θv≤±35°時，均方根誤差小于0.01。此外，在擬合曲線中，一次項a1和常數(shù)項a0與各組數(shù)據(jù)的θv均值具有顯著的相關(guān)性，如圖9所示。隨著傳感器逐漸偏離垂直觀測，一次項a1逐漸減小，與3.2.1 節(jié)TOA 反射率與θv的關(guān)系相反，且a1=b1× cos(θv+b2)+b3；但常數(shù)項a0隨著傳感器逐漸偏離垂直觀測，其數(shù)值逐漸增大，與3.2.1 節(jié)TOA 反射率與θv的關(guān)系相同，同時有：a0=c1× cos(θv+c2)+c3。因此，可以得出TOA反射率與θs及θv的關(guān)系式：

圖9 一次項a1和常數(shù)項a0與觀測天頂角的關(guān)系Fig.9 The relationship between the first order term a1 and the constant term a0 and view zenith angle

式中，θv為觀測天頂角，函數(shù)f1(θv)及f2(θv)描述了TOA 反射率隨θv的變化規(guī)律，cos(θs)是采用控制變量法推導(dǎo)出的TOA反射率隨θs的變化規(guī)律。

3.2.3 散射角對TOA反射率的影響

為刻畫“太陽—目標—傳感器”觀測幾何中方位角差異帶來的TOA 反射率變化，本文引入了Ross?Li 模型體散射核（Zhou 等，2020）中定義的散射角并分析其對TOA 反射率的影響規(guī)律?；讦萻、φs、θv及φv計算散射角Θ的表達式為

式中，θs為太陽天頂角；φs為太陽方位角；θv為傳感器觀測天頂角；φv為傳感器觀測方位角；Θ是散射角，即太陽入射方向與傳感器觀測方向在4π空間中的夾角（Liou，2002）。

針對格爾木沙地穩(wěn)定目標，當(dāng)?35°<θv?θs<35°時，Θ<52°。此時，AQUA/MODIS 傳感器太陽反射率波段的TOA 反射率隨Θ的變化如圖10 所示。其中，后向散射（相對方位角RAA <180°）時Θ記為正，前向散射（相對方位角RAA >180°）時Θ記為負。從圖中可以看出，Θ絕對值越大，太陽入射方向與傳感器觀測方向間的夾角越大，TOA反射率越大。受瑞利散射影響，Θ的變化對藍波段（B3）TOA 反射率的影響最顯著。本文為建TOA 反射率與Θ的關(guān)系模型，使用式（4）所示的傅里葉級數(shù)對不同波段的TOA 反射率數(shù)據(jù)進行擬合，本文選擇n=1。

圖10 AQUA/MODIS傳感器過境格爾木沙地時TOA反射率隨散射角角的變化規(guī)律Fig.10 The TOA reflectance changes with scattering angle when AQUA/MODIS sensor passes through the Golmud Desert

3.3 大氣季節(jié)性影響建模

本文從ECMWF 再分析資料中提取了AQUA/MODIS 傳感器過境格爾木沙地穩(wěn)定目標時對應(yīng)的大氣參數(shù)，其中水汽及臭氧含量變化如圖11 所示，有明顯的季節(jié)性特征。2010 年—2020 年11 年間晴空條件下格爾木沙地實驗區(qū)水汽柱含量呈“M型”分布，春季和秋季水汽柱含量最大，最高為2.05 g/cm2；夏季和冬季水汽柱含量明顯降低，最低為0.04 g/cm2。臭氧柱含量剛好相反，呈“W 型”分布，春季和秋季臭氧柱含量最小，最低為5.11×10?4g/cm2；夏季和冬季臭氧柱含量明顯升高，最高為8.41×10?4g/cm2。2010 年—2020 年格爾木沙地研究區(qū)晴空條件下總AOD@550 nm 無季節(jié)變化特征，但是黑炭、有機質(zhì)、海鹽、鹽酸鹽4 種AOD的變化呈“W 型”的季節(jié)循環(huán)特征，如圖12所示。11 年間，這4 種AOD 從春季到夏季迅速增加，秋季降低，冬季又開始上升。

圖11 2010年—2020年格爾木沙地的大氣季節(jié)性變化規(guī)律Fig.11 Atmospheric seasonal variation from 2010 to 2020 in Golmud Desert

圖12 2010年—2020年格爾木沙地的AOD季節(jié)性變化規(guī)律Fig.12 Seasonal variation of AOD from 2010 to 2020 in Golmud Desert

大氣中臭氧、水汽、氣溶膠光學(xué)厚度等參數(shù)季節(jié)性變化的同時將對TOA 反射率數(shù)據(jù)產(chǎn)生季節(jié)性影響。如圖13 所示，消除θs、θv及Θ等觀測幾何對TOA 反射率的影響后，模型殘差表現(xiàn)出季節(jié)性變化。其中，模型殘差用觀測值與模擬值的比值表示。本文通過使用正弦函數(shù)對不同波段的模型殘差值進行擬合，建立TOA 反射率受臭氧含量、水汽含量、氣溶膠光學(xué)厚度等參數(shù)影響的大氣模型。擬合結(jié)果發(fā)現(xiàn)，盡管B1、B2、B3、B4、B5、B7 波段的振幅不同，但w均為0.0173，即周期約363.19天，近似1年。

圖13 大氣季節(jié)性變化與TOA反射率的關(guān)系Fig.13 The relationship between seasonal variation of the atmosphere and TOA reflectance

3.4 TOA反射率模型

基于3.2 及3.3 節(jié)的分析，建立的TOA 反射率模型如式（5）所示。

4 結(jié)果及驗證

4.1 基于AQUA/MODIS數(shù)據(jù)的建模結(jié)果

由于本文所采用的ERA5 數(shù)據(jù)庫和EAC4 數(shù)據(jù)庫空間分辨率分別為0.25°×0.25°、0.75°×0.75°，因此盡管在數(shù)據(jù)下載及檢查階段已經(jīng)剔除掉15 km×15 km 范圍內(nèi)云覆蓋率較高、過境時刻降水或降雪的影像，但在0.25°× 0.25°的大區(qū)域范圍內(nèi)仍有云覆蓋、降水及降雪。受昆侖山常年積雪影響，格爾木沙地實驗區(qū)0.25°× 0.25°分辨率的夏季雪密度大于0.1 g/cm3。故本文在建立式（5）所示的TOA反射率模型時除了制定與觀測幾何相關(guān)的模型適用性約束條件外，還基于ECMWF 再分析資料制定了大氣參數(shù)等適用性約束條件以篩選MODIS TOA反射率，如表1所示。

表1 與觀測幾何及大氣參數(shù)相關(guān)的TOA反射率模型適用性約束條件Table 1 Observation geometry and atmospheric parameter constraints

本文按照表1制定的觀測幾何及約束條件，篩選2010 年—2019 年AQUAMODIS 的多角度TOA 反射率數(shù)據(jù)共288 景。對建立的TOA 反射率模型，采用最小二乘法，以模擬值與觀測值的誤差最小為約束條件，求解TOA 反射率模型系數(shù)。圖14 展示了B7 波段的模型系數(shù)并對比了該模型模擬TOA反射率與AQUA/MODIS 的B7 波段的觀測TOA 反射率。從圖14 中可以看出：TOA 反射率模型不僅能有效計算不同觀測幾何、地表大氣條件的TOA 反射率，同時保留了陸表穩(wěn)定目標TOA 反射率的季節(jié)性變化特征。二者的相對誤差頻數(shù)分布直方圖如圖15 所示，總體優(yōu)于4%，基本符合正態(tài)分布。統(tǒng)計所有波段對應(yīng)相對誤差的均值、標準差及均方根誤差，結(jié)果見表2。各波段的平均相對誤差在?0.10%以內(nèi)，模型計算值與衛(wèi)星觀測值之間不存在明顯的系統(tǒng)偏差。相對誤差的標準差在3.07%以內(nèi)，均方根誤差小于0.0084。

表2 模型模擬TOA反射率與AQUA/MODIS觀測TOA反射率的平均差異Table 2 The average difference between TOA reflectance simulated by model and observed by AQUA/MODIS

圖14 模型模擬TOA反射率與衛(wèi)星觀測TOA反射率Fig.14 TOA reflectance simulated by model and observed by satellite

圖15 相對誤差頻數(shù)分布直方圖Fig.15 Histogram of relative error frequency distribution

4.2 基于Sentinel-2 及Landsat 8 觀測數(shù)據(jù)的模型驗證

4.2.1 光譜匹配因子

利 用Sentinel?2B/MSI 及Landsat 8 OLI觀測數(shù)據(jù)對本文建立的TOA 反射率模型進行驗證，為解決傳感器間的光譜差異，本文使用式（6）及式（7）將預(yù)測的MODIS傳感器相應(yīng)波段TOA反射率轉(zhuǎn)換成目標傳感器對應(yīng)通道的TOA反射率。

4.2.2 驗證結(jié)果

2018年1月至2021年8月Sentinel?2A/B衛(wèi)星過境格爾木的晴空影像共79 景，過境時1.21°<θv<1.28°、42°<φv<49°、19°<θs<60°、123°<φs<166°，其中30 景（Sentinel?2A 15 景，Sentinel?2B 15 景）符合3.4 節(jié)定義的觀測幾何及大氣參數(shù)約束條件。比較格爾木沙地穩(wěn)定目標模型計算的TOA 反射率（B8A 波段）與Sentinel?2A/B MSI 衛(wèi)星觀測TOA 反射率，如圖16 所示。該波段模型模擬TOA 反射率與衛(wèi)星觀測TOA 反射率間的相對差異如圖17 所示，基本在±2%以內(nèi)，其平均相對誤差為?0.38%，RMSE 約0.0039，在Sentinel?2A/B MSI 傳感器3%的定標不確定度范圍內(nèi)。同時由于地表/大氣的季節(jié)性變化，相對誤差存在1.22%的隨機變化。

圖16 Sentinel?2A/B MSI傳感器B8A模擬TOA反射率與觀測TOA反射率對比Fig.16 The simulated TOA reflectances are compared with observed TOA reflectance of Sentinel?2A/B MSI in B8A

圖17 Sentinel?2A/B MSI傳感器B8A模擬TOA反射率與觀測TOA反射率的相對差異Fig.17 The relative error between simulated TOA reflectances and observed TOA reflectance of Sentinel?2A/B MSI in B8A

圖18 和圖19 展示了本文TOA 反射率模型應(yīng)用于Landsat 8 OLI 傳感器B5 波段的結(jié)果，通過分析可知，相對誤差基本在±2%以內(nèi)，模型計算值與衛(wèi)星觀測值間的RMSE 為0.0047，平均相對誤差為0.71%，相對誤差的標準差為1.35%，表明模型計算值與Landsat 8 OLI 傳感器觀測值具有較好的一致性。

圖18 Landsat 8 OLI傳感器B5波段模擬TOA反射率與觀測TOA反射率對比Fig.18 The simulated TOA reflectances are compared with observed TOA reflectances of Landsat 8 OLI in B5

圖19 Landsat 8 OLI傳感器B5波段模擬TOA反射率與觀測TOA反射率的相對差異Fig.19 The relative error between simulated TOA reflectances and observed TOA reflectances of Landsat 8 OLI in B5

本文提出的TOA 反射率模型計算的TOA 反射率與衛(wèi)星觀測TOA反射率間的相對差異見表3。從表中可知，Sentinel?2B/MSI傳感器在B3、B4、B8A 波段的觀測TOA 反射率與模擬TOA 反射率的平均相對誤差分別為：?0.41%、?0.04%、?0.38%，RMSE 在0.0039 以內(nèi)，相對誤差的標準差不超過1.42%。Landsat 8 OLI 傳感器在B3、B4、B5 波段的觀測TOA 反射率與模擬TOA 反射率的平均相對誤差分別為：1.77%、0.21%、0.71%，RMSE 在0.0055 以內(nèi)，相對誤差的標準差不超過1.72%。由于MSI、OLI 與MODIS 傳感器在短波紅外波段顯著的光譜響應(yīng)函數(shù)差異，導(dǎo)致TOA 反射率模型在Sentinel2A/B MSI、Landsat 8 OLI傳感器B12、B7波段的模擬TOA 反射率與觀測TOA 反射率的RMSE 最高，分別為0.0080、0.0086。對于藍色波段，2 個傳感器的模擬TOA 反射率與觀測TOA 反射率的平均相對差異分別為1.44%、1.50%。

表3 模型預(yù)測值與衛(wèi)星載荷觀測值間相對差異Table 3 The mean and STD of relative error and RMSE between simulated TOA reflectance and observed TOA reflectance

5 結(jié)論

針對國產(chǎn)衛(wèi)星高精度、高頻次在軌輻射定標需求，本文提出了一種綜合考慮“太陽—目標—傳感器”觀測幾何及場地地表/大氣季節(jié)性變化的穩(wěn)定目標TOA反射率模型構(gòu)建方法。以格爾木沙地穩(wěn)定目標為例，結(jié)合Aqua/MODIS 數(shù)據(jù)和ECMWF 的大氣參數(shù)產(chǎn)品，建立了陸表穩(wěn)定目標TOA反射率模型，并以具有高定標精度的Sentinel?2A/B MSI 和Landsat 8 OLI傳感器為參考，對其太陽反射率波段進行了絕對輻射定標驗證。研究結(jié)果表明：本文建立的陸表穩(wěn)定目標（格爾木沙地）TOA 反射率模型較為穩(wěn)定，充分考慮了中國本土陸表穩(wěn)定目標的季節(jié)性特征，能夠提供晴空條件下、滿足觀測幾何及大氣參數(shù)約束條件的光學(xué)衛(wèi)星傳感器不同通道的TOA 反射率參考。其中，模擬的參考TOA 反射率與Sentinel?2A/B MSI 觀測值的平均相對誤差優(yōu)于1.44%，RMSE 在0.0080 以內(nèi)，相對誤差的標準差不超過1.59%；與Landsat 8 OLI 觀測值的平均相對誤差優(yōu)于1.77%，RMSE 在0.0086 以內(nèi)，相對誤差的標準差不超過2.11%。

本文提出的陸表穩(wěn)定目標TOA 反射率模型不僅能用于衛(wèi)星傳感器場地絕對輻射定標，還可用于基于高精度參考衛(wèi)星載荷的交叉輻射定標，這兩種定標方法無需任何同步測量，對國產(chǎn)陸地衛(wèi)星載荷的高精度、高頻次輻射定標具有重要意義。在下一步工作中，本研究將繼續(xù)解決以下問題：（1）分析TOA反射率模型的不確定度以及基于該模型的絕對輻射定標和交叉定標的不確定度。（2）將模型推廣至中國其他陸表穩(wěn)定目標。

志 謝文中使用的MODIS、Sentinel 及Landsat 8 數(shù)據(jù)分別來源于NASA 的LAADS DAAC網(wǎng)站、ESA 網(wǎng)站、USGS 數(shù)據(jù)分發(fā)網(wǎng)站，文中使用的ECMWF 再分析數(shù)據(jù)來源于Atmosphere Data Store 網(wǎng)站及Climate Data Store 網(wǎng)站，在此表示感謝！