用劈窗算法反演地表溫度的通道問題討論

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

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

用劈窗算法反演地表溫度的通道問題討論

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

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

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

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

0 引言

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

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

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

1.1 Planck 函數

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

式中: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 為熱力學溫度，K。

1.2 熱紅外輻射傳輸方程

Chandrasekhar［8］把傳感器接收的熱紅外輻射表達為地物輻射、大氣的上行輻射和大氣下行輻射經地表反射輻射3項之和。熱紅外輻射傳輸方程的建立主要從大氣透過情況、大氣上下行輻射、觀測時衛星天頂角、視場內目標的復雜度和探測器的通道響應等多方面進行考慮。

2 劈窗算法的推導

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

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

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

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

若存在

則式(5)變為

令

則有

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

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

3MODTRAN4模擬和IASI實測數據對比

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

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

從圖1可以看出，MODTRAN模擬和IASI實測的譜線在形態和值的大小上都非常接近，說明用MODTRAN進行相關模擬分析的結果是比較可信的。

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

4.1 大氣貢獻輻射亮度的近似計算

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

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

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

4.2 通道區間Planck函數的線性分析

利用MODTRAN 4模擬了在中緯度夏季大氣模式下，假定海洋表面發射率為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 個通道區間內，波數和輻射亮度間具有較好的線性關系，這種線性關系不會隨著溫度、大氣柱狀水汽含量的改變而改變。

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

在劈窗算法的推導過程中，若 ?B(v，Ts)/?v，/?v 和 ?I(v，T)/?v 存在近似相等關系，則有助于輻射傳輸方程的簡化。觀察這3項發現，其主要差異在于溫度。為了分析?B(v，T)/?v隨溫度的變化趨勢，分別計算其在溫度275 K，280 K，285 K，290 K，295 K，300 K，305 K，310 K 和 315 K時的值。以275 K時的計算結果為基點，分別求與其他溫度狀態下的計算結果差值，即 ?B280，…，315(v，T)/?v－ ?B275(v，T)/?v。以 ?B(v，T)/?v 的差值 δ作為縱坐標軸變量，以波數作為橫坐標軸變量，不同溫度下的?B(v，T)/?v變化趨勢如圖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)區間的通道因 ?B(v，Ts)/?v?v和 ?I(v，T)/?v的近似相等導致的誤差會小于位于8 ～9.09 μm(1250 ～1100 cm－1)區間的通道。

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

為了分析大氣輻射的影響，我們在假定下墊面溫度接近絕對零度、大氣模式選擇中緯度夏季條件下，利用MODTRAN4模擬大氣的分光輻射亮度，如圖5所示。

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

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

4.5 γ值的求解

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

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

圖6 γ的迭代計算Fig.6 γ iterative calculation

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

對于 A組，溫度為 295 K、波數為 909 cm－1(11 μm)時的Planck函數值為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

對于 B組，溫度為295 K、波數為1111 cm－1(9 μm)時的Planck函數值為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

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

5 結論

針對劈窗算法基本公式推導過程中所提出的問題，結合星載探測器的通道響應函數和大氣輻射傳輸模式的模擬數據進行了討論，得出如下結論:

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

2)在所研究的波長范圍內，相鄰通道的大氣輻射亮度非常接近，利用這一特點可以最大可能地削弱大氣影響;

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

［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，等.用陸地衛星TM6數據演算地表溫度的單窗算法［J］.地理學報，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

中國科學院上海技術物理研究所創新專項(編號:Q－ZY－19)及國家科技支撐計劃課題(編號:2012BAH31B02)共同資助。

10.6046/gtzyyg.2012.04.03

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

(責任編輯:刁淑娟)