孟 鵬，胡 勇，鞏彩蘭，栗 琳

(1.中國(guó)科學(xué)院上海技術(shù)物理研究所，上海 200083;2.中國(guó)科學(xué)院研究生院，北京 100049)

用劈窗算法反演地表溫度的通道問(wèn)題討論

孟 鵬1，2，胡 勇1，鞏彩蘭1，栗 琳1，2

(1.中國(guó)科學(xué)院上海技術(shù)物理研究所，上海 200083;2.中國(guó)科學(xué)院研究生院，北京 100049)

劈窗算法是基于熱紅外遙感數(shù)據(jù)反演地表溫度中應(yīng)用較為廣泛、且簡(jiǎn)單有效的算法之一，所使用的熱紅外通道主要位于10～13.3 μm(1000～750 cm－1)波長(zhǎng)區(qū)間內(nèi)，很少考慮8 ～9.09 μm(1250 ～1100 cm－1)區(qū)間內(nèi)的通道數(shù)據(jù)。為了探討更多的適合于反演地表溫度的通道數(shù)據(jù)，結(jié)合劈窗算法基本公式的推導(dǎo)過(guò)程，歸納出了與通道設(shè)置相關(guān)的問(wèn)題，并針對(duì)這些問(wèn)題在10 ～13.3 μm(1000 ～750 cm－1)和8 ～9.09 μm(1250 ～1100 cm－1)區(qū)間內(nèi)進(jìn)行了數(shù)值模擬分析。結(jié)果顯示，用基于此推導(dǎo)劈窗算法，通過(guò)迭代求解，利用8～9.09 μm(1250～1100 cm－1)和10～13.3 μm(1000～750 cm－1)波長(zhǎng)區(qū)間內(nèi)數(shù)據(jù)進(jìn)行溫度反演的結(jié)果非常接近，因此認(rèn)為，將 8 ～9.09 μm(1250 ～1100 cm－1)波長(zhǎng)區(qū)間內(nèi)的數(shù)據(jù)用于劈窗算法反演地表溫度具有一定的可行性。

劈窗算法;地表溫度;熱紅外通道;MODTRAN

0 引言

地球表面溫度是一項(xiàng)重要的地球物理參數(shù)。人們希望獲得不同空間分辨率、任意覆蓋和多時(shí)相的地表溫度數(shù)據(jù)，為此，遙感應(yīng)用工作者提出了一系列用熱紅外遙感數(shù)據(jù)反演地表溫度的算法，如:針對(duì)只有1個(gè)熱紅外通道數(shù)據(jù)(TM/ETM+6)的單通道算法［1，2］、針對(duì)具有 2 個(gè)熱紅外通道數(shù)據(jù)(NOAA/AVHRR)的劈窗算法［3］、針對(duì)具有5個(gè)熱紅外通道數(shù)據(jù)(Terra/ASTER)的溫度發(fā)射率分離算法(TES)［4］以及針對(duì) 7個(gè)紅外通道數(shù)據(jù)(Terra/MODIS)的晝夜算法［5］，等等。單通道算法是基于大氣校正的溫度反演方法，計(jì)算復(fù)雜且誤差較大;TES算法和MODIS晝/夜算法雖然也具有一定的優(yōu)勢(shì)，但需要多個(gè)熱紅外通道數(shù)據(jù)，且計(jì)算復(fù)雜度較高，因而其應(yīng)用受到一定的限制。與單通道和多通道算法相比，針對(duì)2個(gè)熱紅外通道數(shù)據(jù)的劈窗算法計(jì)算簡(jiǎn)單、方法成熟，是目前應(yīng)用最為廣泛的溫度反演算法之一。

當(dāng)前劈窗算法所使用的通道主要位于10～13.3 μm(1000 ～750 cm－1)波長(zhǎng)區(qū)間內(nèi)，很少涉及8～9.09 μm(1250 ～1100 cm－1)區(qū)間內(nèi)的熱紅外通道。從遙感探測(cè)器的設(shè)計(jì)角度考慮，研制用于10 μm內(nèi)長(zhǎng)波紅外探測(cè)的長(zhǎng)線列或面陣器件的技術(shù)成熟度要高，所以在此基礎(chǔ)上進(jìn)行10 μm內(nèi)的熱紅外通道遙感應(yīng)用分析就顯得比較有意義。本文的研究目標(biāo)就是對(duì)8～9.09 μm(1250 ～1100 cm－1)區(qū)間內(nèi)的熱紅外通道用于劈窗算法進(jìn)行地球表面溫度反演的可行性進(jìn)行分析。

1 熱紅外輻射傳輸基本理論

1.1 Planck 函數(shù)

1900年，Planck用量子物理的新概念補(bǔ)充了經(jīng)典物理理論，給出了可以準(zhǔn)確表達(dá)黑體熱紅外輻射能量的Planck函數(shù)［6］。以波數(shù)v(波長(zhǎng)的倒數(shù))表示黑體的分光輻射亮度［7］，其表達(dá)式為

式中:B(v，T)代表黑體的分光輻射亮度，W·cm－2·sr－1·cm;h為普朗克常量，h=6.6252×10－34J·s;k 是玻爾茲曼常量，k=1.38044 ×10－23J·K－1;c是光速，c=2.997925 ×108m·s－1;T 為熱力學(xué)溫度，K。

1.2 熱紅外輻射傳輸方程

Chandrasekhar［8］把傳感器接收的熱紅外輻射表達(dá)為地物輻射、大氣的上行輻射和大氣下行輻射經(jīng)地表反射輻射3項(xiàng)之和。熱紅外輻射傳輸方程的建立主要從大氣透過(guò)情況、大氣上下行輻射、觀測(cè)時(shí)衛(wèi)星天頂角、視場(chǎng)內(nèi)目標(biāo)的復(fù)雜度和探測(cè)器的通道響應(yīng)等多方面進(jìn)行考慮。

2 劈窗算法的推導(dǎo)

為了探討將8 ～9.09 μm(1250 ～1100 cm－1)區(qū)間內(nèi)的數(shù)據(jù)用于劈窗算法進(jìn)行地球表面溫度反演的可行性，同時(shí)降低分析的難度和保證準(zhǔn)確度，選擇發(fā)射率均一且輻射特性接近黑體的海洋表面作為研究對(duì)象，忽略大氣下行輻射的影響。

Mcminllin［9］假定海洋表面發(fā)射率為 1，天空晴好，并且在忽略多次散射的條件下，給出了熱紅外輻射傳輸方程表達(dá)式為

式中:I(v，T)表示大氣頂層輻射亮度;v表示波數(shù);T表示亮度溫度;B是普朗克分光輻射亮度;Ts是地球表面溫度;τ是透過(guò)率;T(P)表示大氣壓強(qiáng)為P時(shí)的大氣溫度;P0是地球表面壓強(qiáng);θ是衛(wèi)星天頂角。根據(jù)積分中值定理，有

式中:I(v，T)，B(v，Ts)和均為含普朗克函數(shù)項(xiàng)。假定在一定的波譜范圍內(nèi)，可以對(duì)普朗克函數(shù)項(xiàng)進(jìn)行泰勒級(jí)數(shù)一階展開，整理后得到

若存在

則式(5)變?yōu)?/p>

令

則有

式中:I1(vr，T1)，I2(vr，T2)分別為 2 個(gè)通道的光譜輻射亮度;γ為和大氣狀態(tài)緊密相關(guān)的量。若γ值確定，則可以確定B(v，Ts)，從而計(jì)算出Ts。

要使上述推導(dǎo)過(guò)程中所作的假設(shè)成立，需要符合5個(gè)條件:①通道所在波長(zhǎng)區(qū)間的普朗克函數(shù)要滿足積分中值定理的要求;②通道所在波段的普朗克函數(shù)具有較好的線性度，可以進(jìn)行泰勒級(jí)數(shù)一階展開;③所選通道對(duì)應(yīng)的 ?B(v，Ts)/?v(v，)/?v和 ?I(v，T)/?v的值近似相等; ④通道內(nèi)的大氣平均輻射亮度Ba(v)近似相等;⑤具備有效的方法求解γ值。

3MODTRAN4模擬和IASI實(shí)測(cè)數(shù)據(jù)對(duì)比

IASI(infrared atmospheric sounding interferometer，紅外大氣探測(cè)干涉儀)是搭載在歐洲極軌氣象衛(wèi)星METOP－A上的探測(cè)器，有8461個(gè)通道，以0.25 cm－1的光譜分辨率覆蓋了3.7 ～15.5 μm 波長(zhǎng)區(qū)間，是目前可以獲得的比較好的紅外光譜數(shù)據(jù)源［10］。以獲取IASI數(shù)據(jù)時(shí)的大氣參數(shù)作為MODTRAN模擬時(shí)的輸入?yún)?shù)，將模擬結(jié)果與IASI實(shí)測(cè)數(shù)據(jù)相比較，如圖1所示。

圖1 MODTRAN 4模擬(左)和IASI實(shí)測(cè)(右)譜線對(duì)比Fig.1 Simulation of MODTRAN(left)vs.Measurements of IASI(right)

從圖1可以看出，MODTRAN模擬和IASI實(shí)測(cè)的譜線在形態(tài)和值的大小上都非常接近，說(shuō)明用MODTRAN進(jìn)行相關(guān)模擬分析的結(jié)果是比較可信的。

4 劈窗算法的通道特性模擬分析

4.1 大氣貢獻(xiàn)輻射亮度的近似計(jì)算

大氣輻射傳輸方程中的積分項(xiàng)表示大氣貢獻(xiàn)的輻射能量。計(jì)算通道區(qū)間在2 μm～1 cm、溫度為233.5 K，273.5 K，310 K 和 373.5 K 時(shí)的 Planck 函數(shù)曲線如圖2所示。

圖2 不同溫度下的Planck函數(shù)曲線Fig.2 Planck curves with different temperatures

從圖2可以看出，Planck函數(shù)曲線滿足中值定理的條件，可以用中值定理近似表達(dá)積分項(xiàng)所代表的大氣貢獻(xiàn)的輻射亮度。

4.2 通道區(qū)間Planck函數(shù)的線性分析

利用MODTRAN 4模擬了在中緯度夏季大氣模式下，假定海洋表面發(fā)射率為1、海洋表面溫度為300 K，大氣柱狀水汽含量為5.0 g/cm2條件下的大氣頂層輻射亮度曲線，如圖3所示。

圖3 通道光譜輻射亮度的線性分析Fig.3 Linear analysis of channel spectral radiance

從圖 3可以看出，在 10～13.3 μm(1000～750 cm－1)和 8 ～9.09 μm(1250 ～1100 cm－1)2 個(gè)通道區(qū)間內(nèi)，波數(shù)和輻射亮度間具有較好的線性關(guān)系，這種線性關(guān)系不會(huì)隨著溫度、大氣柱狀水汽含量的改變而改變。

4.3 ?B(v，T)/?v隨溫度的變化

在劈窗算法的推導(dǎo)過(guò)程中，若 ?B(v，Ts)/?v，/?v 和 ?I(v，T)/?v 存在近似相等關(guān)系，則有助于輻射傳輸方程的簡(jiǎn)化。觀察這3項(xiàng)發(fā)現(xiàn)，其主要差異在于溫度。為了分析?B(v，T)/?v隨溫度的變化趨勢(shì)，分別計(jì)算其在溫度275 K，280 K，285 K，290 K，295 K，300 K，305 K，310 K 和 315 K時(shí)的值。以275 K時(shí)的計(jì)算結(jié)果為基點(diǎn)，分別求與其他溫度狀態(tài)下的計(jì)算結(jié)果差值，即 ?B280，…，315(v，T)/?v－ ?B275(v，T)/?v。以 ?B(v，T)/?v 的差值 δ作為縱坐標(biāo)軸變量，以波數(shù)作為橫坐標(biāo)軸變量，不同溫度下的?B(v，T)/?v變化趨勢(shì)如圖4所示。

圖4 ?B(v，T)/?v隨溫度的變化Fig.4 Changes of ?B(v，T)/?v with different temperatures

從圖4可以看出，靠近800 cm－1處的?B(v，T)/?v隨溫度的變化比較小，位于10 ～13.3 μm(1000 ～750 cm－1)區(qū)間的通道因 ?B(v，Ts)/?v?v和 ?I(v，T)/?v的近似相等導(dǎo)致的誤差會(huì)小于位于8 ～9.09 μm(1250 ～1100 cm－1)區(qū)間的通道。

4.4 通道大氣輻射亮度(v)

為了分析大氣輻射的影響，我們?cè)诩俣ㄏ聣|面溫度接近絕對(duì)零度、大氣模式選擇中緯度夏季條件下，利用MODTRAN4模擬大氣的分光輻射亮度，如圖5所示。

圖5 通道大氣輻射亮度Fig.5 Channel atmospheric radiance

從圖5可以看出，在熱紅外大氣窗區(qū)內(nèi)的10～13.3 μm(1000 ～750 cm－1)和 8 ～9.09 μm(1225 ～1100 cm－1)波長(zhǎng)內(nèi)，大氣輻射亮度隨波數(shù)的變化起伏較小，即在2個(gè)通道設(shè)置比較靠近的條件下，圖中虛線表示的是ASTER 5個(gè)熱紅外通道(B10—B14)的通道響應(yīng)函數(shù)。相鄰?fù)ǖ篱g的大氣輻射亮度的差是比較小的，從而可以近似得到的關(guān)系。

4.5 γ值的求解

Mcmillin提出的利用迭代法求解γ的基本思想是:已知2個(gè)通道或者1個(gè)通道2個(gè)角度的實(shí)測(cè)輻射亮度I1和I2，同時(shí)給定一個(gè)目標(biāo)處的輻射亮度Bsi，然后通過(guò)熱紅外輻射傳輸模式得到2個(gè)計(jì)算輻射亮度I1i和I2i(i表示迭代次數(shù));將Bsi，I1i和 I2i代入式(11)，計(jì)算 γi; 再將 γi，I1和 I2代入式(11)，計(jì)算新的Bsi;重復(fù)上述過(guò)程，直至Bsi的值前(Bsi_a)后(Bsi_b)之差小于某一個(gè)閾值，迭代過(guò)程結(jié)束(圖6)。

為了驗(yàn)證上述迭代法求解γ的效果，以MODTRAN作為大氣輻射傳輸模擬工具，利用ASTER熱紅外探測(cè)器5個(gè)通道的響應(yīng)函數(shù)來(lái)分析迭代求值模型。選 擇 ASTER B11，B12，B13 和 B14，將B 13和B14作為A組，以909 cm－1(11 μm)作為

圖6 γ的迭代計(jì)算Fig.6 γ iterative calculation

參考波數(shù);將B11和B12作為B組，以1111 cm－1(9 μm)作為參考波數(shù)。選擇模擬時(shí)的大氣模式為中緯度夏季大氣(MLS)，下墊面溫度為295 K，其他參數(shù)設(shè)為默認(rèn)值。

對(duì)于 A組，溫度為 295 K、波數(shù)為 909 cm－1(11 μm)時(shí)的Planck函數(shù)值為1.074772E-5 W·cm－2·sr－1·cm;B13的光譜輻射亮度I2=9.63986E-6 W·cm－2·sr－1·cm;B14的光譜輻射亮度I1=1.045429E-5 W·cm－2·sr－1·cm，如表1所示。

表1 ASTER B13和B14的γ迭代①Tab.1 γ iteration of ASTER B13 and B14

對(duì)于 B組，溫度為295 K、波數(shù)為1111 cm－1(9 μm)時(shí)的Planck函數(shù)值為7.271688E-6W·cm－2·sr－1·cm;B11的光譜輻射亮度I2=5.734018E-6 W·cm－2·sr－1·cm;B12的光譜輻射亮度I1=6.688973E-6 W·cm－2·sr－1·cm，如表2所示。

表2 ASTER B11和B12的γ迭代①Tab.2 γ iteration of ASTER B11 and B12

通過(guò)對(duì)A組和B組的數(shù)值模擬分析，盡管2組的參考波數(shù)不同，但都在較少的迭代次數(shù)下收斂。其中，有2個(gè)比較關(guān)鍵的問(wèn)題:①如何給出迭代初始值Bsi;②計(jì)算I1i和I2i所使用輻射傳輸模式。這二者關(guān)系到迭代時(shí)收斂的速度。在上述模擬時(shí)，將I1作為迭代初始值，同時(shí)使用MODTRAN 4作為計(jì)算I1i和I2i的輻射傳輸模式。

5 結(jié)論

針對(duì)劈窗算法基本公式推導(dǎo)過(guò)程中所提出的問(wèn)題，結(jié)合星載探測(cè)器的通道響應(yīng)函數(shù)和大氣輻射傳輸模式的模擬數(shù)據(jù)進(jìn)行了討論，得出如下結(jié)論:

1)以Planck函數(shù)為核函數(shù)的大氣輻射積分項(xiàng)滿足積分中值定理的條件，可以進(jìn)行Planck函數(shù)近似計(jì)算，且在 10 ～13.3 μm(1000 ～750 cm－1)和8 ～9.09 μm(1250 ～1100 cm－1)區(qū)間內(nèi)，Planck 函數(shù)與波數(shù)之間存在較好的線性關(guān)系;

2)在所研究的波長(zhǎng)范圍內(nèi)，相鄰?fù)ǖ赖拇髿廨椛淞炼确浅＝咏?，利用這一特點(diǎn)可以最大可能地削弱大氣影響;

3) 雖然 ?B(v，Ts)/?v和 ?I(v，T)/?v在10 ～13.3 μm(1000 ～750 cm－1)區(qū)間內(nèi)的近似相等關(guān)系優(yōu)于8～9.09μm(1250～1100 cm－1)區(qū)間的，但是利用迭代法求解γ值和輻射亮度對(duì)此并不敏感，所要考慮的是迭代初始值和計(jì)算效率問(wèn)題。因此，將8 ～9.09 μm(1250 ～1100 cm－1)波長(zhǎng)區(qū)間通道獲取的熱紅外遙感數(shù)據(jù)用于劈窗算法進(jìn)行地表溫度反演具有一定的可行性。

［1］Schott J R，Volchok W J.Thematic Mapper Thermal Infrared Calibration［J］.Photogrametric Engineering and Remote Sensing，1985，51(9):1351 －1357.

［2］覃志豪，Zhang M H，Karnieli A，等.用陸地衛(wèi)星TM6數(shù)據(jù)演算地表溫度的單窗算法［J］.地理學(xué)報(bào)，2001，56(4):456 －466.Qin Z H，Zhang M H，Karnieli A，et al.Mono－ window Algorithm for Retrieving Land Surface Temperature from Landsat TM6 Data［J］.Acta Geographica Sinica，2001，56(4):456 － 466(in Chinese with English Abstract).

［3］Price J C.Land Surface Temperature Meas urements from the Split Window Channels of the NOAA7 Advanced Very High Resolution Radiometer［J］.Journal of Gephysical Research，1984，89(D5):7231－7237.

［4］Gillespie A，Rokugawa S，Matsunaga T，el al.Temperature and Emissivity Seperation Algorithm for Advanced Spaceborne Thermal Emission and Reflection Radiometer(ASTER)Images［J］.IEEE Transactions on Geoscience and Remote Sensing，1998，36(4):1113－1126.

［5］Wan Z M，Li Z L.A Physics－ based Algorithm for Retrieving Land Surface Emissivity and Temperature from EOS/MODIS Data［J］.IEEE Transactions on Geoscience and Remote Sensing，1997，35(4):980－996.

［6］Planck M.On the Law of Distribution of Energy in the Normal Spectrum［J］.Annalender Physik，1901，4:553.

［7］Michels T E.Planck Functions and Integrals:Methods of Computation［R］.Washingtion D C:National Aeronautics and Space Administration，1968.

［8］Chandrasekhar S.Radiative Transfer［M］.New York:Dover Publication，1960.

［9］Mcmillin L.Estimation of Sea Surface Temperatures from Two Infrared Window M － easures with Different Absorption［J］.Joural of Geophysical Research，1975，80(C36):5113 －5117.

［10］CNES and EUMETSAT.IASI Program［EB/OL］.http://smsc.cnes.fr/IASI/.

Discussions on Using Channels of Split－window Algorithm to Retrieve Earth Surface Temperature

MENG Peng1，2，HU Yong1，GONG Cai－ lan1，LI Lin1，2
(1.Shanghai Institute of Technical Physics，CAS，Shanghai 200083，China;2.Graduate School of Chinese Academy of Sciences，Beijing 100049，China)

Being simple and effective，the split－window algorithm based on thermal infrared remote sensing is widely used to retrieve surface temperature.The method mainly uses thermal infrared bands in 10 ～ 13.3 μm(1000 ～750 cm－1)range，neglecting bands in 8 ～9.09 μm(1250 ～1100 cm－1)range.This paper analyzes the process of deriving the formula of the split－window algorithm，summarizes the problems associated with the channel setting and makes numerical simulation analysis in the 10 ～13.3 μm(1000 ～750 cm－1)and 8 ～9.09 μm(1250 ～1100 cm－1)ranges to solve the problems.The results show that split－window algorithm derived on the basis of this approach has similar performance in both 10 ～13.3 μm(1000 ～750 cm－1)and 8 ～9.09 μm(1250 ～1100 cm－1)spectral ranges.Therefore，it can be concluded that the spectral range in 8 ～ 9.09 μm(1250 ～ 1100 cm－1)range can also be used to derive split－ window algorithm for thermal remote sensing.

split－window algorithm;Earth surface temperature;thermal infrared band;MODTRAN

TP 79

A

1001－070X(2012)04－0016－05

2011－12－19;

2012－03－08

中國(guó)科學(xué)院上海技術(shù)物理研究所創(chuàng)新專項(xiàng)(編號(hào):Q－ZY－19)及國(guó)家科技支撐計(jì)劃課題(編號(hào):2012BAH31B02)共同資助。

10.6046/gtzyyg.2012.04.03

孟 鵬(1983－)，碩士研究生，主要從事紅外遙感數(shù)據(jù)應(yīng)用研究。E－mail:bushhouse@126.com。

(責(zé)任編輯:刁淑娟)