盧新玉，謝國(guó)輝，李 楊，陳蜀江，馮志敏

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

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

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

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

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

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

0 引言

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

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

1 反演方法

1．1 分裂窗算法

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

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

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

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

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

1．2．1 亮度溫度的計(jì)算

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

MODIS第31、32波段的輻射強(qiáng)度值計(jì)算公式為

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

計(jì)算得到圖像的熱輻射強(qiáng)度之后，便可用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)波段的有效中心波長(zhǎng)，由于第31和第 32波段的波長(zhǎng)區(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)波段的熱輻射強(qiáng)度，由式(5)給出。

1．2．2 大氣透過率的計(jì)算

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

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

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

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

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

2 算法實(shí)現(xiàn)

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

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

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

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

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

2．2 定標(biāo)處理和亮溫計(jì)算

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

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

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

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

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

2．4 云檢測(cè)處理

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

2．5 積雪面積提取

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

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

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

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

2．7 幾何糾正

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

3 反演實(shí)例

以新疆北疆為研究對(duì)象，所用數(shù)據(jù)為新疆氣象局氣候中心提供的MODIS 1B數(shù)據(jù)。對(duì)比驗(yàn)證的雪面溫度觀測(cè)數(shù)據(jù)來自新疆氣象局信息中心。使用設(shè)計(jì)完成的雪面溫度反演程序，以2008年1月18日13時(shí)過境的TERRA數(shù)據(jù)進(jìn)行反演計(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ù)重復(fù)的bow－tie現(xiàn)象;

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

(3)進(jìn)行云檢測(cè);

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

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

4 結(jié)果分析與精度評(píng)價(jià)

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

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

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

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

5 結(jié)論與討論

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

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

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

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

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

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

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

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

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

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

［10］趙英時(shí)，陳冬梅，李小文，等．遙感應(yīng)用分析原理與方法［M］．北京:科學(xué)出版社，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］．北京:機(jī)械工業(yè)出版社，2006:143－152．

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

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

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

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

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

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

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－)，男，碩士研究生，主要研究方向?yàn)橘Y源與環(huán)境遙感監(jiān)測(cè)。

TP 79

A

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

2010－03－13;

2010－04－21

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