基于MODIS/AQUA 數(shù)據(jù)的地表短波凈輻射反演

應(yīng)王敏 ，吳 驊 ※

（1. 中國(guó)科學(xué)院地理科學(xué)與資源研究所/資源與環(huán)境信息系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100101；2. 中國(guó)科學(xué)院大學(xué)/資源與環(huán)境學(xué)院，北京100049）

0 引言

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

基于統(tǒng)計(jì)經(jīng)驗(yàn)關(guān)系的方法是指不考慮太陽(yáng)、大氣與陸地之間復(fù)雜大氣福射傳輸物理機(jī)制過(guò)程，直接建立地面短波輻射的測(cè)量值與遙感觀測(cè)值的經(jīng)驗(yàn)關(guān)系或統(tǒng)計(jì)關(guān)系，以此來(lái)估算地面的短波輻射物理量［3-4］。Pink 和Corio［5］通過(guò)研究發(fā)現(xiàn)，NOAA5 衛(wèi)星的大氣頂層觀測(cè)數(shù)據(jù)與地表凈輻射值具有很強(qiáng)的相關(guān)性，這種關(guān)系隨后也被許多學(xué)者所證實(shí)［6-7］。經(jīng)驗(yàn)統(tǒng)計(jì)法具有步驟簡(jiǎn)單的優(yōu)點(diǎn)，但由于其受到復(fù)雜大氣狀況和復(fù)雜地表特征的限制，普適性較差。

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

基于物理意義的參數(shù)化方法則模擬了大氣與地表間輻射傳輸?shù)奈锢頇C(jī)制，包括大氣反射、散射和吸收作用、地表大氣交互作用等［11］，這些相互作用在地表短波輻射反演中進(jìn)行的參數(shù)化構(gòu)建［12］。例如，Pinker 和Ewing［13］對(duì)大氣層中水汽、云、氣溶膠的吸收與散射作用進(jìn)行參數(shù)化處理，成功反演得到地表短波凈輻射，并分析出太陽(yáng)天頂角（SZA）對(duì)地表短波凈輻射影響較大?；谖锢硪饬x的參數(shù)化方法具有一定的物理意義，能夠較為精確地估算模擬結(jié)果，適用性較強(qiáng)。但其模型構(gòu)建過(guò)程相對(duì)復(fù)雜，模型精度受大氣地表輸入?yún)?shù)的精確度和獲取困難度的影響。此外，大部分參數(shù)化方案是基于比較理想的前提假設(shè)下提出的，部分參數(shù)的取值具有區(qū)域性和經(jīng)驗(yàn)性。

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

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

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

1.1 研究區(qū)概況

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

圖1 SURFRAD 地面實(shí)測(cè)站點(diǎn)分布Fig.1 The distribution of SURFRAD ground stations

1.2 MODIS/AQUA 數(shù)據(jù)

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

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

1.3 MODTRAN 模擬數(shù)據(jù)

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

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

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

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

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

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

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

由于現(xiàn)實(shí)情況下，我們往往無(wú)法獲取到大氣頂層向上輻射通量（Fu）這個(gè)物理量，考慮到Fu與大氣頂層寬通道反照率（r）有轉(zhuǎn)換關(guān)系：

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

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

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

3 地表短波凈輻射反演模型構(gòu)建與驗(yàn)證

由上可知，構(gòu)建地表短波輻射反演模型需要兩個(gè)步驟：（1）窄通道表觀反射率向?qū)捦ǖ婪凑章兽D(zhuǎn)換；（2）利用大氣頂層向上輻射通量與地表短波凈輻射的非線性關(guān)系反演NSSR。公式（5）～（7）是第一個(gè)步驟所需要擬合的公式，其擬合系數(shù)與太陽(yáng)天頂角SZA有關(guān)，每一個(gè)VZA 對(duì)應(yīng)一套系數(shù)（28 個(gè)系數(shù)）。這里共模擬了36 個(gè)SZA 值，故我們的系數(shù)查找表的維度為36*32，由于篇幅關(guān)系，這里并不作展示。

圖2 展示了步驟一的擬合結(jié)果，橫坐標(biāo)表示MODTRAN 模擬出來(lái)的寬波段反射率，縱坐標(biāo)是我們基于窄波段反射率模擬得到的寬波段反射率?？梢园l(fā)現(xiàn)，散點(diǎn)的樣本數(shù)為13.996 8 萬(wàn)個(gè)，偏差Bias 為0，均方根誤差RMSE 為0.011，決定系數(shù)R2為0.997，該結(jié)果保證了大氣頂層窄波段反射率到大氣頂層寬波段反照率轉(zhuǎn)換的方法精度。圖中散點(diǎn)幾乎平均分布在1∶1 參考線附近，并且每類代表性地物具有特定的反照率范圍。

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

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

表2 NSSR 反演模型構(gòu)建步驟二——方程（1-3）系數(shù)擬合統(tǒng)計(jì)Table 2 Fitting statistics of equation coefficients

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

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

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

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

表3 地表短波凈輻射在SURFRAD 站點(diǎn)的反演誤差統(tǒng)計(jì)Table 3 Error statistics of NSSR retrieval at seven SURFRAD sites

4 結(jié)論

地表短波凈輻射是地表輻射平衡中的一個(gè)重要組分，是許多全球、區(qū)域氣候、水文、陸面過(guò)程模型中的重要參數(shù)，其高精度遙感反演對(duì)全球和區(qū)域氣候變化、能量平衡、生態(tài)環(huán)境、地面模型構(gòu)建、大氣循環(huán)探究等領(lǐng)域具有相當(dāng)重要的實(shí)際意義。該文以MODIS/AQUA 數(shù)據(jù)為驅(qū)動(dòng)源，構(gòu)建了具有一定泛化能力的高精度地表短波凈輻射反演模型。

在反演模型構(gòu)建過(guò)程中，利用大氣傳輸模型MODTRAN 模擬了13.996 8 萬(wàn)個(gè)種復(fù)雜情況下的地表短波輻射?；诖罅康哪M數(shù)據(jù)集，擬合地表短波凈輻射與大氣頂層觀測(cè)值的非線性關(guān)系，該反演模型分為兩個(gè)步驟：（1）窄通道表觀反射率向?qū)捦ǖ婪凑章兽D(zhuǎn)換；（2）利用大氣頂層向上輻射通量與地表短波凈輻射的非線性關(guān)系反演NSSR。步驟一構(gòu)建結(jié)果表明，其偏差Bias 為0，均方根誤差RMSE 為0.011，決定系數(shù)R2為0.997。步驟二的偏差Bias 為0 W/m2，均方根誤差RMSE 為7.29 W/m2，決定系數(shù)R2為0.999。誤差統(tǒng)計(jì)結(jié)果說(shuō)明地表短波凈輻射反演模型的構(gòu)建精度較高。

構(gòu)建完地表短波凈輻射反演模型后，利用地表輻射觀測(cè)網(wǎng)（SURFRAD）7 個(gè)地面實(shí)測(cè)站點(diǎn)數(shù)據(jù)，結(jié)合MODIS/AQUA 遙感數(shù)據(jù)，對(duì)反演模型進(jìn)行驗(yàn)證。驗(yàn)證結(jié)果表明，其總體偏差Bias 為-4.8 W/m2，均方根誤差RMSE 為77.1 W/m2，決定系數(shù)R2為0.9106。研究發(fā)現(xiàn)，超70%樣本的反演誤差低于50 W/m2，在不同站點(diǎn)具有不同的性能表現(xiàn)。