基于MODIS/AQUA 數據的地表短波凈輻射反演

應王敏 ，吳 驊 ※

（1. 中國科學院地理科學與資源研究所/資源與環境信息系統國家重點實驗室，北京100101；2. 中國科學院大學/資源與環境學院，北京100049）

0 引言

地表短波凈輻射（Net Surface Shortwave Radiation，NSSR），指地表吸收的短波輻射通量（0.3～5 μm）。地表短波凈輻射是地表輻射平衡中的一個重要組分，是氣候形成及氣候變化的主要依據，影響著物質輸送、地氣能量交換、水汽流動、生物光合作用等［1-2］。根據地表短波凈輻射遙感反演方法的不同，可以分為基于統計經驗關系的方法、基于物理過程的方法、基于物理意義的參數化方法以及混合模型方法。

基于統計經驗關系的方法是指不考慮太陽、大氣與陸地之間復雜大氣福射傳輸物理機制過程，直接建立地面短波輻射的測量值與遙感觀測值的經驗關系或統計關系，以此來估算地面的短波輻射物理量［3-4］。Pink 和Corio［5］通過研究發現，NOAA5 衛星的大氣頂層觀測數據與地表凈輻射值具有很強的相關性，這種關系隨后也被許多學者所證實［6-7］。經驗統計法具有步驟簡單的優點，但由于其受到復雜大氣狀況和復雜地表特征的限制，普適性較差。

基于物理過程的方法指根據大氣輻射傳輸理論，給定大氣和地表參數，通過求解復雜的輻射傳輸方程計算得到地表短波輻射［8］。根據太陽輻射經大氣到達地表的傳輸過程，Joseph 等［9］提出了一種快速但精確的計算大氣吸收—散射輻射通量的方法，Delta-Eddington 近似法，通過聯合Dirac delta 函數及二流近似，解決了Eddington 近似處理高度不對稱相函數精度不高的問題。Lenoble［10］討論了不同的大氣多次散射輻射傳輸方法，還比較模型的復雜性及精度。該方法的優點是具有明確的物理意義，計算精確。但其由于輸入參數眾多，計算復雜，難以直接應用于像元級運算。

基于物理意義的參數化方法則模擬了大氣與地表間輻射傳輸的物理機制，包括大氣反射、散射和吸收作用、地表大氣交互作用等［11］，這些相互作用在地表短波輻射反演中進行的參數化構建［12］。例如，Pinker 和Ewing［13］對大氣層中水汽、云、氣溶膠的吸收與散射作用進行參數化處理，成功反演得到地表短波凈輻射，并分析出太陽天頂角（SZA）對地表短波凈輻射影響較大。基于物理意義的參數化方法具有一定的物理意義，能夠較為精確地估算模擬結果，適用性較強。但其模型構建過程相對復雜，模型精度受大氣地表輸入參數的精確度和獲取困難度的影響。此外，大部分參數化方案是基于比較理想的前提假設下提出的，部分參數的取值具有區域性和經驗性。

混合模型法是將經驗統計法和物理模型參數化法結合使用，既考慮了大氣輻射傳輸方程的物理機制，又發揮了經驗統計方法的間接操作和較高精度的優勢。例如Tang 等［14］利用輻射傳輸模型MODTRAN 模擬生成各種復雜情況下的短波凈輻射，通過統計回歸的方式構建大氣頂層觀測數據與地表短波凈輻射的非線性關系，以此來實現MODIS/TERRA 地表短波凈輻射的遙感反演。該類方法具有一定物理意義，反演精度高，模型構建相對簡潔等優勢，是目前反演地表短波凈輻射廣泛應用的方法之一［15-17］。

文章基于混合模型法，試圖在大量MODTRAN 模擬集上構建針對MODIS/AQUA 數據的地表短波凈輻射反演模型，并利用具有不同下墊面情況的多站點長時間序列數據對反演模型進行綜合全面的魯棒性分析和精度檢驗，最終發展出具有一定泛化能力的適用于MODIS/AQUA 數據的高精度地表短波凈輻射反演模型。

1 研究區與數據

1.1 研究區概況

地表輻射觀測網（SURFRAD）是一個常用的地表輻射觀測網絡，由美國國家海洋和大氣局提供（http：//www.srrb.noaa.gov/surfrad）。該地表輻射觀測網絡現有7 個觀測站點，其空間分布情況如圖1 所示，分別位于蒙大拿州（Fort Peck，MT；北緯48.31°，西經105.10°）、卡羅拉多州（Table Mountain，CO；北緯40.13°，西經105.24°）、伊利羅伊州（Bondville，IL；北緯40.05°，西經88.37°）、密西西比州（Goodwin Creek，MS；北緯34.25°，西經 89.87°）、賓夕法尼亞（Penn State，PA；北緯 40.72°，西經 77.93°）、內華達州（Desert Rock，NV；北緯36.62°，西經116.02°）和南達科他州（Sioux Falls，SD；北緯43.73°，西經96.62°）。站點位置選擇考慮了美國的不同氣候區特點、下墊面類型的異質性，相應的測量值能夠有效地代表該區域的特征。SURFRAD 觀測網絡為美國本土的氣候、數值預報及衛星遙感監測等研究提供準確、連續、長期的地面輻射觀測數據，包括地表短波下行輻射、地表短波上行輻射、地表短波凈輻射等，其采樣的時間分辨率為1 min（2009 年1 月1 日之前的數據記錄間隔為3 min）。此外，SURFRAD 還提供了質量控制（QC）、大氣溫度、相對濕度、風速、風向、氣壓等輔助變量。

圖1 SURFRAD 地面實測站點分布Fig.1 The distribution of SURFRAD ground stations

1.2 MODIS/AQUA 數據

搭載于美國AQUA 太陽同步極軌衛星平臺上的中分辨率成像光譜儀（MODIS），擁有36 個通道，覆蓋了可見光—紅外整個波譜［18］。MODIS/AQUA 在地方時下午13∶30 過境，提供了相應MYD 系列產品，空間分辨率一般為1 km。NASA 研究人員在原始輻亮度數據產品的基礎上生成了很多種不同級別、不同類型的遙感應用產品（陸地產品、大氣產品、海洋產品等），提供具有時空連續性的地表監測資料，為地表短波輻射的研究提供有效的手段，滿足地球科學多方面的研究需求［19］。

該文所用到的產品包括MYD02，MYD03，MYD05 和MYD35，下載了過境研究區域7 個地面站點的2017 年全年遙感影像（https：//ladsweb.modaps.eosdis.nasa.gov/）。其中MYD02 產品可以提供通道輻射亮度信息，MYD03 產品提供了觀測幾何參數，MYD05 產品可提供大氣可降水量變量，而MYD35 產品則提供相應的云掩膜信息，可用于判斷晴空和有云條件。由于受儀器自身硬件的影響，MYD02 的第6 個通道存在部分數據缺失的現象，在大氣層的遙感影像上表現出較多的條紋，因此該文僅利用了MYD02 的前5 個通道和第7 個通道來反演地表的短波凈輻射。

1.3 MODTRAN 模擬數據

大氣輻射傳輸模型MODTRAN 5 可以有效模擬各種典型地表、大氣剖面、氣溶膠、云、觀測幾何等復雜條件下的地表短波輻射。基于MODTRAN 的模擬就能獲取不同地表和大氣狀況下的模擬數據，從而構建地表短波輻射反演模型。

設置模擬的波譜范圍為300～5 000 nm，光譜分辨率為2 波數（cm-1）。為了充分模擬各種復雜情況，從MODTRAN 默認數據庫里選擇了9 種有代表性的地類波譜反射率數據作為輸入，代表著不同下墊面類型，包括草地、濕地、土壤、森林、沙漠、城市、海水、新雪、海冰。在模擬時，采用了MODTRAN 5 提供的6 種標準大氣廓線（熱帶大氣，中緯度夏季大氣，中緯度冬季大氣，副極帶夏季大氣，副極帶冬季大氣以及美國1976 年標準大氣US76）。同時，使用了3 種氣溶膠類型（鄉村氣溶膠、海洋性氣溶膠以及對流層氣溶膠），相應的能見度設置為缺省值（能見度取值范圍為10～35 km）。除了晴空天氣以外，還考慮了3 種云模型，即層云、高層云和積云。考慮到一般MODIS 數據的最大觀測天頂角小于 65°，故模擬了 6 種觀測天頂角 VZA（0°，33.5°，44.4°，51.3°，56.2°，60°）。此外，考慮了36 種太陽天頂角SZA 變化，即0°～70°，間隔2°。將這些設置參數條件排列組合，總共模擬了139968 個樣本量的模擬數據集。

MODTRAN 模型可以模擬相應情況下的波譜輻射亮度（存儲在TP7 文件里）、多層高度對應的向上輻射通量和向下輻射通量（存儲在FLX 文件里）等。地表短波凈輻射（NSSR）可以通過將FLX 文件中最底層高度對應的向下輻射通量減去向上輻射通量得到；MODIS 通道輻射亮度（Li）可以將TP7 文件中波譜輻射亮度與MODIS 通道響應函數卷積得到。

2 地表短波凈輻射反演方法

利用大氣頂層TOA 向上輻射通量（Fu）與地表短波凈輻射（NSSR）的非線性關系可有效反演NSSR［14］，該方法已經被廣泛應用到各種衛星傳感器，并得到了較好的反演結果。公式（1）是該方法的核心思想，NSSR 與Fu整體上呈現一種負相關的關系：

式（1）中，μs表示太陽天頂角SZA 的余弦值，D是地球與太陽的天文距離（與日期有關），E0是進入大氣頂部的太陽輻照度（可以通過MODTRAN 內部文件獲取，取1 368 W/m2）。此外，α′和β′是方程系數，可以通過公式（2）和公式（3）獲得：

式（2）～（3）中，μs表示太陽天頂角SZA 的余弦值，w代表著大氣可降水含量，a1～a7，x，y，z是方程待擬合系數，與天氣情況有關（晴空或有云）。聯合公式（1）、（2）、（3），基于上述139968 對NSSR 和Fu的樣本集，結合MATLAB 和最小二乘法，可以擬合出特定天氣情況下的適用系數a1～a7，x，y，z。

由于現實情況下，我們往往無法獲取到大氣頂層向上輻射通量（Fu）這個物理量，考慮到Fu與大氣頂層寬通道反照率（r）有轉換關系：

式（4）中，μs，D，E0的物理意義和公式（1）一致。因此若能獲取到大氣頂層寬通道反照率（r），則可以將其轉換成大氣頂層向上輻射通量（Fu），進而得到NSSR。但是，大部分衛星傳感器只提供了窄波段的表觀反射率，所以需將窄通道表觀反射率向寬通道反照率進行轉換，通過將傳感器幾個窄波段反射率進行線性組合進而獲取寬波段反照率。以MODIS/AQUA 為例，寬波段反照率可以表示為：

式（5）～（7）中，i表示波段，ρi表示該波段的表觀反射率，Li表示該波段的輻射亮度；μv表示衛星觀測天頂角VZA 的余弦值；μs表示太陽天頂角SZA 的余弦值；D是地球與太陽的天文距離；方程（6）中的c1i，c2i，c3i，c4i是系數，需要通過最小二乘法進行擬合；Ei是該波段大氣頂層的平均太陽輻照度，可以通過MODTRAN 內部文件和MODIS/AQUA 通道響應函數卷積得到，表1 列出了MODIS/AQUA 選定的6 個通道響應的平均太陽輻照度。

表1 MODIS/AQUA 選定的6 個通道大氣頂層太陽輻照度Table 1 the TOA solar irradiance of the six selected MODIS/AQUA radiance bands

3 地表短波凈輻射反演模型構建與驗證

由上可知，構建地表短波輻射反演模型需要兩個步驟：（1）窄通道表觀反射率向寬通道反照率轉換；（2）利用大氣頂層向上輻射通量與地表短波凈輻射的非線性關系反演NSSR。公式（5）～（7）是第一個步驟所需要擬合的公式，其擬合系數與太陽天頂角SZA有關，每一個VZA 對應一套系數（28 個系數）。這里共模擬了36 個SZA 值，故我們的系數查找表的維度為36*32，由于篇幅關系，這里并不作展示。

圖2 展示了步驟一的擬合結果，橫坐標表示MODTRAN 模擬出來的寬波段反射率，縱坐標是我們基于窄波段反射率模擬得到的寬波段反射率。可以發現，散點的樣本數為13.996 8 萬個，偏差Bias 為0，均方根誤差RMSE 為0.011，決定系數R2為0.997，該結果保證了大氣頂層窄波段反射率到大氣頂層寬波段反照率轉換的方法精度。圖中散點幾乎平均分布在1∶1 參考線附近，并且每類代表性地物具有特定的反照率范圍。

圖2 NSSR 反演模型構建步驟一——大氣頂層寬波段反照率與估算值的散點Fig.2 Comparison of the estimated TOA broadband albedos with MODTRAN albedos

反演模型的第二個步驟，是先利用公式（4）將大氣頂層寬波段反照率轉換到大氣頂層向上輻射通量，再通過公式（1）～（3）擬合地表短波凈輻射與大氣頂層向上輻射通量的非線性關系。步驟二中的方程系數與天氣狀況（晴空或有云）有關，即晴空或有云各自有一套方程系數，表2 是通過最小二乘法和MATLAB 擬合得到的系數結果。圖3展示了NSSR 反演模型步驟二的結果，橫坐標表示MODTRAN 模擬出來的地表短波凈輻射，縱坐標是基于大氣頂層向上輻射通量得到的地表短波凈輻射。晴空和有云散點的總樣本數為13.996 8 萬個，偏差Bias 為0 W/m2，均方根誤差RMSE 為7.29 W/m2，決定系數R2為0.999。晴空狀況的NSSR 數值范圍為0～1 200 W/m2，而有云情況的數值范圍為0～400 W/m2，可以發現有云狀況時反演結果相對比較離散。由此可見，基于MODTRAN模擬數據的地表短波凈輻射反演模型的構建精度較高。

表2 NSSR 反演模型構建步驟二——方程（1-3）系數擬合統計Table 2 Fitting statistics of equation coefficients

圖3 NSSR 反演模型構建二——地表短波凈輻射和估算值的散點Fig.3 Comparison of the estimated NSSR with MODTRAN-modeled NSSR

基于MODTRAN 模擬數據構建出的地表短波輻射反演模型，需要經過地面站點實測數據來檢驗。進行地面站點驗證前，為了減少異常值對數據驗證結果的影響，將地面站點部分存在異常情況的數據進行了剔除。異常值的判斷原則是，某瞬時測量值的前后5 min 內的標準差大于給定閾值，這樣可以消除現場測量時隨機噪聲等帶來的不良影響。以2017 年全年為研究時間，需要提取MODIS/AQUA 產品中距離相應地面站點最近的遙感影像像元值。提取的變量信息包括：MYD02 產品中的通道輻射亮度Li，MYD03 產品中的經緯度、太陽天頂角SZA、觀測天頂角VZA 等觀測幾何參數，MYD05 產品中的大氣可降水量w，MYD35 產品中的云掩膜信息。利用MODIS 云掩膜信息區分晴空和云兩種情況，將Li、SZA、VZA、w等自變量輸入到已構建好的NSSR 反演模型。將反演模型得到的NSSR 和地面實測NSSR 的進行比較，統計反演誤差。

圖4 是NSSR 反演模型的地面驗證結果。圖4（a）展示了站點散點圖，可以看出，模型驗證的總體偏差Bias 為-4.8 W/m2，均方根誤差RMSE 為77.1 W/m2，決定系數R2為0.9106。經過對比探究，發現NSSR 反演模型在不同站點上的性能表現不同，其在“Bondville，IL”站點的驗證結果最好，而在“Sioux Falls，SD”站點則誤差較大（表3）。圖4（b）則是模型誤差分布直方圖，可以發現超70%樣本的反演誤差低于50 W/m2，且誤差符合正態分布。綜上，所構建的地表短波凈輻射反演模型在地面站點驗證結果良好，具有一定泛化能力。

圖4 NSSR 反演模型地面驗證（a）散點圖（b）誤差直方圖Fig.4 Validation of NSSR retrieval model using in situ measurements.

表3 地表短波凈輻射在SURFRAD 站點的反演誤差統計Table 3 Error statistics of NSSR retrieval at seven SURFRAD sites

4 結論

地表短波凈輻射是地表輻射平衡中的一個重要組分，是許多全球、區域氣候、水文、陸面過程模型中的重要參數，其高精度遙感反演對全球和區域氣候變化、能量平衡、生態環境、地面模型構建、大氣循環探究等領域具有相當重要的實際意義。該文以MODIS/AQUA 數據為驅動源，構建了具有一定泛化能力的高精度地表短波凈輻射反演模型。

在反演模型構建過程中，利用大氣傳輸模型MODTRAN 模擬了13.996 8 萬個種復雜情況下的地表短波輻射。基于大量的模擬數據集，擬合地表短波凈輻射與大氣頂層觀測值的非線性關系，該反演模型分為兩個步驟：（1）窄通道表觀反射率向寬通道反照率轉換；（2）利用大氣頂層向上輻射通量與地表短波凈輻射的非線性關系反演NSSR。步驟一構建結果表明，其偏差Bias 為0，均方根誤差RMSE 為0.011，決定系數R2為0.997。步驟二的偏差Bias 為0 W/m2，均方根誤差RMSE 為7.29 W/m2，決定系數R2為0.999。誤差統計結果說明地表短波凈輻射反演模型的構建精度較高。

構建完地表短波凈輻射反演模型后，利用地表輻射觀測網（SURFRAD）7 個地面實測站點數據，結合MODIS/AQUA 遙感數據，對反演模型進行驗證。驗證結果表明，其總體偏差Bias 為-4.8 W/m2，均方根誤差RMSE 為77.1 W/m2，決定系數R2為0.9106。研究發現，超70%樣本的反演誤差低于50 W/m2，在不同站點具有不同的性能表現。