基于IDL的MODIS 1B數(shù)據(jù)雪面溫度反演

盧新玉，謝國輝，李 楊，陳蜀江，馮志敏

(1．新疆氣象局，烏魯木齊 830002;2．中國氣象局烏魯木齊沙漠氣象研究所，烏魯木齊 830002;3．新疆師范大學，烏魯木齊 830054)

基于IDL的MODIS 1B數(shù)據(jù)雪面溫度反演

盧新玉1，2，3，謝國輝1，李 楊2，陳蜀江3，馮志敏1

(1．新疆氣象局，烏魯木齊 830002;2．中國氣象局烏魯木齊沙漠氣象研究所，烏魯木齊 830002;3．新疆師范大學，烏魯木齊 830054)

對MODIS數(shù)據(jù)的分裂窗算法進行了簡要介紹，通過交互式數(shù)據(jù)語言(IDL)編程，實現(xiàn)了直接利用MODIS lB數(shù)據(jù)的雪面溫度反演，并將反演結(jié)果存貯為標準的HDF文件格式，以供其他軟件使用。以我國新疆北部(北疆)為例，將該算法的反演結(jié)果與氣象站雪面溫度觀測資料進行對比，結(jié)果表明:系統(tǒng)反演得到的雪面溫度分布規(guī)律與氣象觀測資料是一致的，反演的平均誤差為1．73℃，基本反映了北疆地區(qū)的雪面溫度分布情況。

分裂窗算法;雪面溫度;反演;MODIS;IDL

0 引言

積雪表面溫度是積雪研究必不可少的參數(shù)之一［1］。它不僅是雪蓋的基本特性，也指示著融雪過程所處的狀態(tài)，甚至可以通過對雪面溫度跟蹤分析融雪徑流的整個過程。積雪表面溫度時空分布對于積雪消融、變質(zhì)及地表大氣能量交換過程均有重要影響，是融雪徑流研究的重要參數(shù)。因此，利用遙感技術(shù)進行雪面溫度反演和動態(tài)監(jiān)測具有重要意義。

分裂窗算法是目前為止發(fā)展最為成熟的地表溫度反演算法，有著廣泛的應(yīng)用［2－4］。本文基于MODIS數(shù)據(jù)的地表溫度分裂窗算法［5］，并結(jié)合積雪特征，利用IDL語言程序設(shè)計，實現(xiàn)雪面溫度的反演算法，為積雪研究過程中快速實時獲取雪面溫度資料提供了便利。

1 反演方法

1．1 分裂窗算法

分裂窗算法主要是針對NOAA/AVHRR的第4、5熱紅外波段提出來的。在MODIS的8個熱紅外波段中，第31、32波段最接近 AVHRR第4、5波段，因而最適用于分裂窗算法［7］。由于熱紅外遙感中的未知變量比波段方程數(shù)多，信息量不足，所以地表溫度反演算法的推導都是根據(jù)一些模擬簡化過程來從波段方程中求解地表溫度變量，不同的簡化方法產(chǎn)生了不同的算法。本文基于覃志豪等［5］提出的適用于MODIS數(shù)據(jù)的地表溫度反演算法，通過IDL編程，實現(xiàn)了直接利用MODIS lB數(shù)據(jù)的雪面溫度反演。該分裂窗算法的表達式為

式中，TS為地表溫度，T31、T32是 MODIS第31、32波段的亮度溫度(單位K)，可用Planck輻射方程獲取;A0、A1和A2是系數(shù)，可用式(2)計算，即

式中，a31、b31、a32和 b32是常量。取 a31=－64．603 63，b31=0．440 817;a32= － 68．725 75，b32=0．473 453，其余中間參數(shù)分別計算如下

式中，τi(θ)是 i(i=31，32)波段視角為 θ時的大氣透過率;εi是i(i=31，32)波段的地表比輻射率。

1．2 分裂窗算法各中間參數(shù)的計算

1．2．1 亮度溫度的計算

MODIS圖像是用DN值表示的，因此，要計算星上亮溫，必須先將熱輻射波段的DN值轉(zhuǎn)換成相應(yīng)的輻射強度值，然后再用Planck函數(shù)求解星上亮溫。

MODIS第31、32波段的輻射強度值計算公式為

式中，Ii是 i(i=31，32)波段的熱輻射強度;DNi是 i(i=31，32)波段的 DN 值;scalei、offseti分別為MODIS第i(i=31，32)波段的輻射定標常量，可從MODIS數(shù)據(jù)集的屬性數(shù)據(jù)中查出。

計算得到圖像的熱輻射強度之后，便可用Planck函數(shù)求解出星上亮度溫度，即

式中，Ti是MODIS第i(i=31，32)波段的亮度溫度，即式(1)中的T31和T32;h為普朗克常數(shù)，取值6．626 ×10－34J·s;c為光速，取值2．998 ×108m/s;λi是第i(i=31，32)波段的有效中心波長，由于第31和第 32波段的波長區(qū)間分別為 10．78～11．28 μm和 11．77 ～12．27 μm，所以，λi可分別取其中值(λ31=11．03 μm 和λ32=12．02 μm);k為玻耳茲曼常數(shù)，取值1．380 6 ×10－23J/K;Ii是 MODIS第i(i=31，32)波段的熱輻射強度，由式(5)給出。

1．2．2 大氣透過率的計算

大氣透過率是分裂窗算法中的一項基本參數(shù)，它主要受大氣水汽含量的影響。對于MODIS數(shù)據(jù)地表溫度反演，大氣水汽含量可通過MODIS第2、19波段來反演，然后再根據(jù)大氣水汽含量與大氣透過率之間的函數(shù)關(guān)系來估計大氣透過率［8］。大氣水汽含量的計算公式為

式中，ω是大氣水汽含量(g/cm2);α、β是常量(α =0．02，β =0．651);ref19、ref2分別是 MODIS 第19、2波段的地面反射率，計算方法見式(9)。

根據(jù)大氣透過率和大氣水汽含量的負線性關(guān)系，毛克彪［6］通過LOWTRAN(大氣模擬軟件)模擬MODIS第31、32波段透過率與大氣水汽含量的關(guān)系，對冬季模擬數(shù)據(jù)進行線性回歸，得到MODIS第31、32波段的大氣透過率與大氣水汽含量間的線性方程分別為

1．2．3 積雪比輻射率確定

地表比輻射率表示地表物質(zhì)的熱特征，其大小主要取決于地表的物質(zhì)結(jié)構(gòu)、組成成分、表面狀態(tài)及物理性質(zhì)等，并隨所測輻射能的波長和觀測角度等條件的變化而變化［9，10］。本文采用 ASTER 光譜庫［11］中提供的積雪光譜確定MODIS第31、32波段積雪比輻射率，分別取0．982 291和0．962 046。

2 算法實現(xiàn)

根據(jù)反演流程(圖1)，采用結(jié)構(gòu)化的方法來設(shè)計程序。將各部分功能均寫為函數(shù)進行調(diào)用，分7個函數(shù)塊實現(xiàn)計算過程。

圖1 雪面溫度反演流程(Bi表示第i波段數(shù)據(jù))Fig．1 The flowchart of retrieving snow surface temperature(Bi means the data of band i)

2．1 MODIS 數(shù)據(jù)讀取

MODIS 1B數(shù)據(jù)采用HDF文件格式，MODIS數(shù)據(jù)通道信息和定位信息都按照相應(yīng)的名稱貯存在科學數(shù)據(jù)集SD(Scientific DataSet)中，讀取這些信息時先要將名稱轉(zhuǎn)換為索引號，然后按照索引號來檢索數(shù)據(jù)。定位信息包括經(jīng)緯度、太陽高度角及方位角等信息。IDL語言中提供了很多函數(shù)用來讀取HDF文件，使得操作HDF文件非常方便［12］。為提高數(shù)據(jù)存貯效率，MODIS 1B產(chǎn)品只將分辨率相同、光譜性質(zhì)相同的數(shù)據(jù)放在同一SDS對象內(nèi)，稱為波段組。讀取通道信息需先清楚所需波段數(shù)據(jù)在哪一SDS對象內(nèi)，具體如表1所示。

表1 MODIS 1B產(chǎn)品中的各波段組Tab．1 The various band groups of MODIS 1B data

下面以第31通道為例，詳細介紹如何利用IDL語言讀取通道數(shù)據(jù)和定標屬性(分號為注釋標記，其后為注釋語句)。

以上程序完成了第31通道數(shù)據(jù)和用于輻射定標計算參數(shù)的讀取。其他通道數(shù)據(jù)、地理數(shù)據(jù)以及傳感器信息數(shù)據(jù)的讀取方式同上。

2．2 定標處理和亮溫計算

將MODIS數(shù)據(jù)的DN值進行定標計算，轉(zhuǎn)換為反射率或輻射度。反射率需要進行太陽高度角校正，即把不同太陽高度角下的探測數(shù)據(jù)轉(zhuǎn)換成相當于天頂時的探測值，校正公式為

式中，R是校正后的反射率;ref是定標后的反射率;θ是太陽高度角。由于MODISl km數(shù)據(jù)中的經(jīng)緯度和太陽高度角數(shù)據(jù)維數(shù)是通道數(shù)據(jù)維數(shù)的1/5，而在同名MOD 03數(shù)據(jù)中的經(jīng)緯度信息和太陽高度角則是和1 km通道數(shù)據(jù)一一對應(yīng)的，因此本文中的太陽高度角數(shù)據(jù)直接從MOD 03數(shù)據(jù)獲取。讀取太陽高度角數(shù)據(jù)的方法為

定標后的熱輻射波段數(shù)據(jù)需要計算亮度溫度，計算方法如1．2．1所述。

2．3 大氣透過率、地表比輻射率的算法實現(xiàn)

根據(jù)前述大氣透過率和地表比輻射率的計算公式設(shè)計函數(shù)，返回大氣透過率、地表比輻射率的二維數(shù)組，在程序中進行調(diào)用。

2．4 云檢測處理

為了提高反演的精度，在進行雪面溫度反演處理前必須對MODIS 1B數(shù)據(jù)進行云檢測，去除那些受云干擾的像元。目前針對MODIS 1B數(shù)據(jù)提出了許多云檢測的方法［13－16］，本文參考陳俊蕙提出的云檢測算法［15］，并根據(jù)積雪下墊面特點對其算法閾值進行了改進。利用云檢測算法最終得到一個二值的數(shù)據(jù)集，0表示云，1表示非云。用此數(shù)據(jù)與積雪溫度數(shù)據(jù)進行乘積運算，得到去除云的積雪溫度圖。

2．5 積雪面積提取

利用歸一化雪蓋指數(shù)算法，根據(jù)文獻［17］所述方法進行積雪面積提取，得到一個二值數(shù)據(jù)集(0表示非積雪，1表示積雪)。將此數(shù)據(jù)與經(jīng)過云檢測的積雪溫度數(shù)據(jù)進行乘積運算，從而得到去除云積雪區(qū)的雪面溫度圖。

2．6 雪面溫度的算法實現(xiàn)與結(jié)果保存

計算出大氣透過率和地表比輻射率這兩個參數(shù)后，積雪溫度就可用分裂窗算法(式(1)～(8))來計算，這在IDL中很容易實現(xiàn)。雪面溫度計算完成后，用IDL中的HDF讀寫函數(shù)，可以把計算結(jié)果保存為HDF格式，具體代碼如下:

云檢測結(jié)果、雪蓋面積和經(jīng)緯度信息都可按上述方法存貯到HDF文件中，分別創(chuàng)建Cloud、SCA、Longitude和Latitude數(shù)據(jù)后，在SnowST之后依次添加。按照該方式存貯的SnowST產(chǎn)品數(shù)據(jù)大大降低了數(shù)據(jù)量，便于數(shù)據(jù)分發(fā)和共享，該格式遵循HDF標準，可以在ENVI等其他商業(yè)軟件中使用。

2．7 幾何糾正

在IDL中調(diào)用ENVI函數(shù)可實現(xiàn)圖像的幾何糾正。主要思路是:讀取MOD 03文件中的Latitude和Longitude數(shù)據(jù)，然后把這個經(jīng)緯度坐標數(shù)組傳遞到ENVI_REGISTER_DOIT函數(shù)的pts參數(shù)，就可以得到幾何糾正的結(jié)果。主要代碼如下:

3 反演實例

以新疆北疆為研究對象，所用數(shù)據(jù)為新疆氣象局氣候中心提供的MODIS 1B數(shù)據(jù)。對比驗證的雪面溫度觀測數(shù)據(jù)來自新疆氣象局信息中心。使用設(shè)計完成的雪面溫度反演程序，以2008年1月18日13時過境的TERRA數(shù)據(jù)進行反演計算，得到北疆地區(qū)雪面溫度分布如圖2所示。

圖2 2008年1月18日北疆地區(qū)雪面溫度分布Fig．2 Snow surface temperature distribution map in north of Xinjiang on January，18，2008

主要步驟如下:

(1)利用modis tools工具下的bow－tie correction命令去除相鄰掃描行之間數(shù)據(jù)重復的bow－tie現(xiàn)象;

(2)進行物理定標，將DN值轉(zhuǎn)化為反射亮度或輻射亮度;

(3)進行云檢測;

(4)用第4、6波段數(shù)據(jù)計算NDSI，提取積雪面積;

(5)幾何糾正。由于MODIS資料提供了精確的經(jīng)緯度信息，所以幾何糾正時完全可以讀取其自帶的經(jīng)緯度數(shù)據(jù)，本文使用IDL程序進行等經(jīng)緯度投影。

4 結(jié)果分析與精度評價

由圖2可以看出，低溫區(qū)主要分布在阿勒泰中部、北部以及伊犁地區(qū);高溫區(qū)主要分布在塔城北部地區(qū)。雪面溫度(注:夏季為草面溫度)觀測項目是近幾年氣象站才開始的一項新內(nèi)容，它是由埋于積雪表面的探頭感應(yīng)測得，然后由自動站采集，是雪面的真實溫度。由于所用衛(wèi)星資料的過境時間為13:14，故本文選用北疆氣象站13:00的觀測數(shù)據(jù)對反演溫度值進行精度驗證，其結(jié)果如表2所示。

表2 北疆地區(qū)實測雪面溫度與反演結(jié)果對比Tab．2 Comparison of retrieved values of snow surface temperature in north of Xinjiang with the observed values

分析表明，所選17個站點的平均誤差為1．73℃，反演誤差在2℃以內(nèi)的占64．7%，在2～3℃的占23．5%，在3℃以上的占11．8%。部分站點誤差相對較大的原因可能是由于圖像像元與地面實測點配準的不確定性造成的。根據(jù)數(shù)理統(tǒng)計原理，對17個站點處的反演溫度和觀測值進行一元線性回歸(圖3)。擬合結(jié)果表明，在95%的置信水平下，兩者之間一元線性相關(guān)性R2=0．901，反演得到的雪面溫度圖基本反映了北疆地區(qū)的雪面溫度分布情況。

圖3 反演值與實測值的一元線性回歸Fig．3 Unary linear regression between retrieved values and measured values

5 結(jié)論與討論

基于分裂窗算法，從業(yè)務(wù)運行的需求出發(fā)，在IDL 7．0環(huán)境下對系統(tǒng)的數(shù)據(jù)處理流程和功能進行了詳細分析，提出了構(gòu)建反演系統(tǒng)的總體構(gòu)架;采用結(jié)構(gòu)化開發(fā)方法對系統(tǒng)框架進行軟件實現(xiàn)，并詳細介紹了IDL語言在實現(xiàn)MODIS lB數(shù)據(jù)雪面溫度反演的編繹方法，最后給出了應(yīng)用實例。結(jié)果表明，本文方法操作簡便、運算速度快，能從海量的MODIS數(shù)據(jù)中直接提取有效數(shù)據(jù)，可用于快速、精確反演地表雪面溫度，從而為積雪的動態(tài)監(jiān)測，以及后續(xù)的融雪徑流、雪水當量研究提供重要依據(jù)。

［1］曹梅盛，李 新，陳賢章，等．冰凍圈遙感［M］．北京:科學出版社，2005:74－80．

［2］武 堅，孟憲紅，呂世華．基于MODIS數(shù)據(jù)的金塔綠洲地表溫度反演［J］．高原氣象，2009，28(3):523 －529．

［3］包 剛，包玉海，李慧靜，等．用MODIS數(shù)據(jù)和分裂窗算法反演內(nèi)蒙古地區(qū)的地表溫度［J］．測繪科學，2009，34(1):32 －34．

［4］許國鵬，李仁東，劉可群，等．基于MODIS數(shù)據(jù)的湖北省地表溫度反演研究［J］．華中師范大學學報:自然科學版，2007，41(1):143－147．

［5］覃志豪，高懋芳，秦曉敏，等．農(nóng)業(yè)旱災監(jiān)測中的地表溫度遙感反演方法——以 MODIS數(shù)據(jù)為例［J］．自然災害學報，2005，14(4):64－71．

［6］孫 亮，孫 睿，賈成剛，等．MODIS數(shù)據(jù)反演地表溫度劈窗算法比較［J］．北京師范大學學報:自然科學版，2008，44(4):434 －438．

［7］毛克彪，覃志豪，施建成，等．針對MODIS影像的劈窗算法研究［J］．武漢大學學報:信息科學版，2005，30(8):703 －707．

［8］毛克彪，覃志豪，王建明，等．針對MODIS數(shù)據(jù)的大氣水汽含量反演及31和32波段透過率計算［J］．國土資源遙感，2005(1):26－30．

［9］覃志豪，李文娟，徐 斌，等．陸地衛(wèi)星TM6波段范圍內(nèi)地表比輻射率的估計［J］．國土資源遙感，2004(3):28－33．

［10］趙英時，陳冬梅，李小文，等．遙感應(yīng)用分析原理與方法［M］．北京:科學出版社，2003:106－110．

［11］ASTER Spectral Library．Reproduced from the ASTER Spectral Library Through the Courtesy of the Jet Propulsion Laboratory［EB/OL］．［2010 －2 －10］．http://speclib．jpl．nasa．gov，2006

［12］閻殿武．IDL可視化工具入門與提高［M］．北京:機械工業(yè)出版社，2006:143－152．

［13］祝善友，尹 球，匡定波．極軌氣象衛(wèi)星圖像的云自動檢測方法研究［J］．遙感信息，2007(1):37 －40．

［14］何全軍，曹 靜，黃 江，等．基于多光譜綜合的MODIS數(shù)據(jù)云檢測研究［J］．國土資源遙感，2006(3):19－22．

［15］陳俊蕙．MODIS數(shù)據(jù)云檢測算法研究［D］．北京:中國農(nóng)業(yè)大學，2007:22－42．

［16］劉玉潔，楊忠東，劉 健，等．MODIS遙感信息處理原理與算法［M］．北京:科學出版社，2001:109－130．

［17］黃 鎮(zhèn)，崔彩霞．基于EOS/MODIS的新疆積雪監(jiān)測［J］．冰川凍土，2006，28(3):343 －347．

(責任編輯:刁淑娟)

The Retrieval of Snow Surface Temperature from MODIS 1B Data Based on IDL

LU Xin － yu1，2，3，XIE Guo － hui1，LI Yang2，CHEN Shu － jiang3，F(xiàn)ENG Zhi－ min1
(1．Xinjiang Meteorological Bureau，Urumqi830002，China;2．Institute of Desert Meteorology，CMA，Urumqi830002，China;3．Xinjiang Normal University，Urumqi830054，China)

In this paper，the splits window algorithm developed by Qin was used to retrieve the snow surface temperature(TS)from MODIS 1B data by using Interactive Data Language(IDL)programming．An IDL program was written for the TS retrieval．The retrieved TS data were saved in the HDF file readable into other software．Details of image processing and programming were shown in this paper．Finally，the authors chose northern Xinjiang as an example to demonstrate the applicability of the approach．The retrieved TS data were compared with the observed snow temperatures from automatic meteorological stations．The result of the comparison indicates that the retrieved TS distribution is consistent with the observed one，and the average error is1．73℃，which is acceptable．Therefore，the TS distribution in northern Xinjiang could be rapidly retrieved using the proposed IDL approach．

Splits window algorithm;Snow surface temperature;Retrieval;MODIS;IDL

盧新玉(1979－)，男，碩士研究生，主要研究方向為資源與環(huán)境遙感監(jiān)測。

TP 79

A

1001－070X(2010)04－0029－05

2010－03－13;

2010－04－21

中央級公益性科研院所基本科研業(yè)務(wù)項目(編號:IDM200704、IDM200808)、沙漠氣象科學研究基金項目(編號:SQJ2009006)及科技部公益性行業(yè)(氣象)科研專項(編號:GYHY200706008)共同資助。