基于IDL的MODIS 1B数据雪面温度反演

卢新玉，谢国辉，李 杨，陈蜀江，冯志敏

(1．新疆气象局，乌鲁木齐 830002;2．中国气象局乌鲁木齐沙漠气象研究所，乌鲁木齐 830002;3．新疆师范大学，乌鲁木齐 830054)

基于IDL的MODIS 1B数据雪面温度反演

卢新玉1，2，3，谢国辉1，李 杨2，陈蜀江3，冯志敏1

(1．新疆气象局，乌鲁木齐 830002;2．中国气象局乌鲁木齐沙漠气象研究所，乌鲁木齐 830002;3．新疆师范大学，乌鲁木齐 830054)

对MODIS数据的分裂窗算法进行了简要介绍，通过交互式数据语言(IDL)编程，实现了直接利用MODIS lB数据的雪面温度反演，并将反演结果存贮为标准的HDF文件格式，以供其他软件使用。以我国新疆北部(北疆)为例，将该算法的反演结果与气象站雪面温度观测资料进行对比，结果表明:系统反演得到的雪面温度分布规律与气象观测资料是一致的，反演的平均误差为1．73℃，基本反映了北疆地区的雪面温度分布情况。

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

0 引言

积雪表面温度是积雪研究必不可少的参数之一［1］。它不仅是雪盖的基本特性，也指示着融雪过程所处的状态，甚至可以通过对雪面温度跟踪分析融雪径流的整个过程。积雪表面温度时空分布对于积雪消融、变质及地表大气能量交换过程均有重要影响，是融雪径流研究的重要参数。因此，利用遥感技术进行雪面温度反演和动态监测具有重要意义。

分裂窗算法是目前为止发展最为成熟的地表温度反演算法，有着广泛的应用［2－4］。本文基于MODIS数据的地表温度分裂窗算法［5］，并结合积雪特征，利用IDL语言程序设计，实现雪面温度的反演算法，为积雪研究过程中快速实时获取雪面温度资料提供了便利。

1 反演方法

1．1 分裂窗算法

分裂窗算法主要是针对NOAA/AVHRR的第4、5热红外波段提出来的。在MODIS的8个热红外波段中，第31、32波段最接近 AVHRR第4、5波段，因而最适用于分裂窗算法［7］。由于热红外遥感中的未知变量比波段方程数多，信息量不足，所以地表温度反演算法的推导都是根据一些模拟简化过程来从波段方程中求解地表温度变量，不同的简化方法产生了不同的算法。本文基于覃志豪等［5］提出的适用于MODIS数据的地表温度反演算法，通过IDL编程，实现了直接利用MODIS lB数据的雪面温度反演。该分裂窗算法的表达式为

式中，TS为地表温度，T31、T32是 MODIS第31、32波段的亮度温度(单位K)，可用Planck辐射方程获取;A0、A1和A2是系数，可用式(2)计算，即

式中，a31、b31、a32和 b32是常量。取 a31=－64．603 63，b31=0．440 817;a32= － 68．725 75，b32=0．473 453，其余中间参数分别计算如下

式中，τi(θ)是 i(i=31，32)波段视角为 θ时的大气透过率;εi是i(i=31，32)波段的地表比辐射率。

1．2 分裂窗算法各中间参数的计算

1．2．1 亮度温度的计算

MODIS图像是用DN值表示的，因此，要计算星上亮温，必须先将热辐射波段的DN值转换成相应的辐射强度值，然后再用Planck函数求解星上亮温。

MODIS第31、32波段的辐射强度值计算公式为

式中，Ii是 i(i=31，32)波段的热辐射强度;DNi是 i(i=31，32)波段的 DN 值;scalei、offseti分别为MODIS第i(i=31，32)波段的辐射定标常量，可从MODIS数据集的属性数据中查出。

计算得到图像的热辐射强度之后，便可用Planck函数求解出星上亮度温度，即

式中，Ti是MODIS第i(i=31，32)波段的亮度温度，即式(1)中的T31和T32;h为普朗克常数，取值6．626 ×10－34J·s;c为光速，取值2．998 ×108m/s;λi是第i(i=31，32)波段的有效中心波长，由于第31和第 32波段的波长区间分别为 10．78～11．28 μm和 11．77 ～12．27 μm，所以，λi可分别取其中值(λ31=11．03 μm 和λ32=12．02 μm);k为玻耳兹曼常数，取值1．380 6 ×10－23J/K;Ii是 MODIS第i(i=31，32)波段的热辐射强度，由式(5)给出。

1．2．2 大气透过率的计算

大气透过率是分裂窗算法中的一项基本参数，它主要受大气水汽含量的影响。对于MODIS数据地表温度反演，大气水汽含量可通过MODIS第2、19波段来反演，然后再根据大气水汽含量与大气透过率之间的函数关系来估计大气透过率［8］。大气水汽含量的计算公式为

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

根据大气透过率和大气水汽含量的负线性关系，毛克彪［6］通过LOWTRAN(大气模拟软件)模拟MODIS第31、32波段透过率与大气水汽含量的关系，对冬季模拟数据进行线性回归，得到MODIS第31、32波段的大气透过率与大气水汽含量间的线性方程分别为

1．2．3 积雪比辐射率确定

地表比辐射率表示地表物质的热特征，其大小主要取决于地表的物质结构、组成成分、表面状态及物理性质等，并随所测辐射能的波长和观测角度等条件的变化而变化［9，10］。本文采用 ASTER 光谱库［11］中提供的积雪光谱确定MODIS第31、32波段积雪比辐射率，分别取0．982 291和0．962 046。

2 算法实现

根据反演流程(图1)，采用结构化的方法来设计程序。将各部分功能均写为函数进行调用，分7个函数块实现计算过程。

图1 雪面温度反演流程(Bi表示第i波段数据)Fig．1 The flowchart of retrieving snow surface temperature(Bi means the data of band i)

2．1 MODIS 数据读取

MODIS 1B数据采用HDF文件格式，MODIS数据通道信息和定位信息都按照相应的名称贮存在科学数据集SD(Scientific DataSet)中，读取这些信息时先要将名称转换为索引号，然后按照索引号来检索数据。定位信息包括经纬度、太阳高度角及方位角等信息。IDL语言中提供了很多函数用来读取HDF文件，使得操作HDF文件非常方便［12］。为提高数据存贮效率，MODIS 1B产品只将分辨率相同、光谱性质相同的数据放在同一SDS对象内，称为波段组。读取通道信息需先清楚所需波段数据在哪一SDS对象内，具体如表1所示。

表1 MODIS 1B产品中的各波段组Tab．1 The various band groups of MODIS 1B data

下面以第31通道为例，详细介绍如何利用IDL语言读取通道数据和定标属性(分号为注释标记，其后为注释语句)。

以上程序完成了第31通道数据和用于辐射定标计算参数的读取。其他通道数据、地理数据以及传感器信息数据的读取方式同上。

2．2 定标处理和亮温计算

将MODIS数据的DN值进行定标计算，转换为反射率或辐射度。反射率需要进行太阳高度角校正，即把不同太阳高度角下的探测数据转换成相当于天顶时的探测值，校正公式为

式中，R是校正后的反射率;ref是定标后的反射率;θ是太阳高度角。由于MODISl km数据中的经纬度和太阳高度角数据维数是通道数据维数的1/5，而在同名MOD 03数据中的经纬度信息和太阳高度角则是和1 km通道数据一一对应的，因此本文中的太阳高度角数据直接从MOD 03数据获取。读取太阳高度角数据的方法为

定标后的热辐射波段数据需要计算亮度温度，计算方法如1．2．1所述。

2．3 大气透过率、地表比辐射率的算法实现

根据前述大气透过率和地表比辐射率的计算公式设计函数，返回大气透过率、地表比辐射率的二维数组，在程序中进行调用。

2．4 云检测处理

为了提高反演的精度，在进行雪面温度反演处理前必须对MODIS 1B数据进行云检测，去除那些受云干扰的像元。目前针对MODIS 1B数据提出了许多云检测的方法［13－16］，本文参考陈俊蕙提出的云检测算法［15］，并根据积雪下垫面特点对其算法阈值进行了改进。利用云检测算法最终得到一个二值的数据集，0表示云，1表示非云。用此数据与积雪温度数据进行乘积运算，得到去除云的积雪温度图。

2．5 积雪面积提取

利用归一化雪盖指数算法，根据文献［17］所述方法进行积雪面积提取，得到一个二值数据集(0表示非积雪，1表示积雪)。将此数据与经过云检测的积雪温度数据进行乘积运算，从而得到去除云积雪区的雪面温度图。

2．6 雪面温度的算法实现与结果保存

计算出大气透过率和地表比辐射率这两个参数后，积雪温度就可用分裂窗算法(式(1)～(8))来计算，这在IDL中很容易实现。雪面温度计算完成后，用IDL中的HDF读写函数，可以把计算结果保存为HDF格式，具体代码如下:

云检测结果、雪盖面积和经纬度信息都可按上述方法存贮到HDF文件中，分别创建Cloud、SCA、Longitude和Latitude数据后，在SnowST之后依次添加。按照该方式存贮的SnowST产品数据大大降低了数据量，便于数据分发和共享，该格式遵循HDF标准，可以在ENVI等其他商业软件中使用。

2．7 几何纠正

在IDL中调用ENVI函数可实现图像的几何纠正。主要思路是:读取MOD 03文件中的Latitude和Longitude数据，然后把这个经纬度坐标数组传递到ENVI_REGISTER_DOIT函数的pts参数，就可以得到几何纠正的结果。主要代码如下:

3 反演实例

以新疆北疆为研究对象，所用数据为新疆气象局气候中心提供的MODIS 1B数据。对比验证的雪面温度观测数据来自新疆气象局信息中心。使用设计完成的雪面温度反演程序，以2008年1月18日13时过境的TERRA数据进行反演计算，得到北疆地区雪面温度分布如图2所示。

图2 2008年1月18日北疆地区雪面温度分布Fig．2 Snow surface temperature distribution map in north of Xinjiang on January，18，2008

主要步骤如下:

(1)利用modis tools工具下的bow－tie correction命令去除相邻扫描行之间数据重复的bow－tie现象;

(2)进行物理定标，将DN值转化为反射亮度或辐射亮度;

(3)进行云检测;

(4)用第4、6波段数据计算NDSI，提取积雪面积;

(5)几何纠正。由于MODIS资料提供了精确的经纬度信息，所以几何纠正时完全可以读取其自带的经纬度数据，本文使用IDL程序进行等经纬度投影。

4 结果分析与精度评价

由图2可以看出，低温区主要分布在阿勒泰中部、北部以及伊犁地区;高温区主要分布在塔城北部地区。雪面温度(注:夏季为草面温度)观测项目是近几年气象站才开始的一项新内容，它是由埋于积雪表面的探头感应测得，然后由自动站采集，是雪面的真实温度。由于所用卫星资料的过境时间为13:14，故本文选用北疆气象站13:00的观测数据对反演温度值进行精度验证，其结果如表2所示。

表2 北疆地区实测雪面温度与反演结果对比Tab．2 Comparison of retrieved values of snow surface temperature in north of Xinjiang with the observed values

分析表明，所选17个站点的平均误差为1．73℃，反演误差在2℃以内的占64．7%，在2～3℃的占23．5%，在3℃以上的占11．8%。部分站点误差相对较大的原因可能是由于图像像元与地面实测点配准的不确定性造成的。根据数理统计原理，对17个站点处的反演温度和观测值进行一元线性回归(图3)。拟合结果表明，在95%的置信水平下，两者之间一元线性相关性R2=0．901，反演得到的雪面温度图基本反映了北疆地区的雪面温度分布情况。

图3 反演值与实测值的一元线性回归Fig．3 Unary linear regression between retrieved values and measured values

5 结论与讨论

基于分裂窗算法，从业务运行的需求出发，在IDL 7．0环境下对系统的数据处理流程和功能进行了详细分析，提出了构建反演系统的总体构架;采用结构化开发方法对系统框架进行软件实现，并详细介绍了IDL语言在实现MODIS lB数据雪面温度反演的编绎方法，最后给出了应用实例。结果表明，本文方法操作简便、运算速度快，能从海量的MODIS数据中直接提取有效数据，可用于快速、精确反演地表雪面温度，从而为积雪的动态监测，以及后续的融雪径流、雪水当量研究提供重要依据。

［1］曹梅盛，李 新，陈贤章，等．冰冻圈遥感［M］．北京:科学出版社，2005:74－80．

［2］武 坚，孟宪红，吕世华．基于MODIS数据的金塔绿洲地表温度反演［J］．高原气象，2009，28(3):523 －529．

［3］包 刚，包玉海，李慧静，等．用MODIS数据和分裂窗算法反演内蒙古地区的地表温度［J］．测绘科学，2009，34(1):32 －34．

［4］许国鹏，李仁东，刘可群，等．基于MODIS数据的湖北省地表温度反演研究［J］．华中师范大学学报:自然科学版，2007，41(1):143－147．

［5］覃志豪，高懋芳，秦晓敏，等．农业旱灾监测中的地表温度遥感反演方法——以 MODIS数据为例［J］．自然灾害学报，2005，14(4):64－71．

［6］孙 亮，孙 睿，贾成刚，等．MODIS数据反演地表温度劈窗算法比较［J］．北京师范大学学报:自然科学版，2008，44(4):434 －438．

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

［8］毛克彪，覃志豪，王建明，等．针对MODIS数据的大气水汽含量反演及31和32波段透过率计算［J］．国土资源遥感，2005(1):26－30．

［9］覃志豪，李文娟，徐 斌，等．陆地卫星TM6波段范围内地表比辐射率的估计［J］．国土资源遥感，2004(3):28－33．

［10］赵英时，陈冬梅，李小文，等．遥感应用分析原理与方法［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］．北京:机械工业出版社，2006:143－152．

［13］祝善友，尹 球，匡定波．极轨气象卫星图像的云自动检测方法研究［J］．遥感信息，2007(1):37 －40．

［14］何全军，曹 静，黄 江，等．基于多光谱综合的MODIS数据云检测研究［J］．国土资源遥感，2006(3):19－22．

［15］陈俊蕙．MODIS数据云检测算法研究［D］．北京:中国农业大学，2007:22－42．

［16］刘玉洁，杨忠东，刘 健，等．MODIS遥感信息处理原理与算法［M］．北京:科学出版社，2001:109－130．

［17］黄 镇，崔彩霞．基于EOS/MODIS的新疆积雪监测［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，FENG 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－)，男，硕士研究生，主要研究方向为资源与环境遥感监测。

TP 79

A

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

2010－03－13;

2010－04－21

中央级公益性科研院所基本科研业务项目(编号:IDM200704、IDM200808)、沙漠气象科学研究基金项目(编号:SQJ2009006)及科技部公益性行业(气象)科研专项(编号:GYHY200706008)共同资助。