MODIS数据地表温度反演分裂窗算法的IDL实现

第29卷第3期2006年06月

测绘与空间地理信息

GEOMATICS&SPATIALINFORMATIONTECHNOLOGY

Vol.29,No.3Jun.,2006

MODIS数据地表温度反演分裂窗

算法的IDL实现

姜立鹏

1,3

,覃志豪

1,2

,谢　雯

3

(1.南京大学国际地球系统科学研究所,江苏南京210093;2.中国农业科学院农业资源与农业区划研究所,

北京100081;3.南京大学城市与资源学系,江苏南京210093)

摘要:地表温度是气象、水文、生态等研究领域中的一个重要参数。本文针对MODIS数据的分裂窗算法进行了

简要的介绍,并对参数的获取进行了分析。义务化流程生产地表温度产品,我们在IDL6.0环境下,,操作简便,不需人为干预就可快速反演地表温度,非常适合批量计算MOIDS关键词:MODIS;地表温度;分裂窗算法;IDL中图分类号:P208　　　文献标识码:B)03-0114-04

AlgorithmtoRetrieveLandSurface

forMODISDataUsingIDL

JIANGLi2peng,QINZhi2hao

1,3

1,2

,XIEWen

3

(1.InternationalInstituteforEarthSystemScienceofNanjingUniversity,Nanjing210093,China;2.InstituteofNaturalResources

andRegionalPlanning,ChineseAcademyofAgriculturalSciences,Beijing100081,China;

3.UrbanandResourceDepartmentofNanjingUniversity,Nanjing210093,China)

Abstract:Landsurfacetemperature(LST)isanimportantparameterinmeteorology,hydrology,ecologyandsoon.Atpresent,MO2DISiswidelyusedforLSTretrieving,andthesplitswindowalgorithmisoneofthemostapprovedmethods.Althoughtheprecisionofthisalgorithmisquitehigh,itisdifficulttooperateonpracticalapplicationbecauseofitscomplexity.Wehaveprogrammedthisalgo2rithmusingIDLlanguage,andrealizedtheLST’sretrievingautomatically.Theresultsindicatethatthisprocedureisfast,easytohan2dle,andcanbeusedforproducingtheLSTproductsinbatches.

Keywords:MODIS;landsurfacetemperature;splitswindowalgorithm;IDL

0　引　言

MODIS(ModerateResolutionImagingSpectroradiome2ter,中等分辨率成像光谱辐射仪)有36个光谱通道,其中29～36通道为热红外通道,可用来监测地球表面的热量

确、简便地实现MODIS数据地表温度反演成为一个迫切需要解决的问题。

Qinetal.(2001)在文献[1]中提出了两因素地表温

度反演分裂窗算法,该算法只需要两个参数并且精度较高。覃志豪等人在以前他本人针对AVHRR提出的劈窗算法基础上,改进并提出了适用于MODIS数据的地表温度反演算法

[2]

变化。

它每天可获得同一地区的白天和夜间的重复观测资料(低纬地区为2d),双星可达到每天4次的遥感数据。对我国白天地表温度反演而言,需要东部、中部和西部3景遥感数据才能覆盖我国全部领土,双星将需要6景遥感数据。而各应用领域(如干旱监测)要求我们必须快速反演出地表温度,以达到各种应用目的。因此,如何快速、精

收稿日期:2005-11-08

基金项目:自然科学基金项目资助(40471096)

,并且算法需要的基本参数都可以从MO2

DIS数据中反演得到。该算法已经被推荐并应用于中国MODIS地表温度产品生产。为了进行义务化流程生产,

减少工作人员的工作量,本文先对这一义务化算法进行介绍,然后对这一算法的程序(IDL)实现进行简要介绍。

作者简介:姜立鹏(1982-),男,硕士生,主要从事热红外遥感理论与应用、数字图像处理、地理信息系统应用研究。

第3期姜立鹏等:MODIS数据地表温度反演分裂窗算法的IDL

实现115

1　分裂窗算法简介及其参数求解

1.1　分裂窗算法

查出;scale32和offset32为MODIS第32波段的辐射定标常量,也可从MODIS数据集的属性数据中查出。

计算得到图像的热辐射强度之后,便可用Planck函数求解出星上亮度温度。对31和32波段分别应用Plank函数,并化简可得31和32波段的亮温计算公式:

(1)

T31=K31,1/ln(1+K31,2/rad31)T32=K32,1/ln(1+K32,2/rad32)

(9)(10)

我们采用的算法是覃志豪在[2]中提出的适用于

MODIS数据的地表温度反演算法,该算法的公式如下:

Ts=A0+A1T31-A2T32

式中,Ts是地表温度(K),T31和T32分别是MODIS第31和

32波段的亮度温度。A0,A1和A2是分裂窗算法的参数,分

式中,K31,1=729.541636,K31,2=1304.413871K,

K32,1=474.684780,K32,2=1196.978785。

别定义如下:

A0=a31D32(1-C31-D31)/(D32C31-D31C32)-a32D31(1-C32-D32)/(D32C31-D31C32)A1=1+D31/(D32C31-D31C32)+

b31D32(1-C31-D31)/(D32C31-D31C32)A2=D31/(D32C31-D31C32)+

b32D31(1-C32-D32)/(D32C31-D3132)

4(3)(2)

1.2.2　大气透过率的计算

大气透过率是分裂窗算法中的一项基本参数,它主要受大气水汽含量的影响。对于数据地表温度反演,I19波段来反演,]

[3]

所示)。大气水汽含量的计

2

:

(11)

在这里,a31,b31,a31,b32是常量,=.=0.440817,a32=-68.,=0.β)w=((α-ln(ref19/ref2))/

-2

式中,w是大气水分含量(g/cm);α和β是常量,分别

τθ)i=εii(

θ)][1+(1-εθ)]Dr=[1-τi(i)τi(

(5)(6)

取α=0.02和β=0.651;ref19和ref2分别是MODIS第

19和2波段的地面反射率,由下式计算:

ref2=scale2(band2-offset2)ref19=scale19(band19-offset19)

(12)(13)

θ)是i(i=31,32)波段视角为θ的大气透过率;其中,τi(

τθ)是i波段视角为θ的大气透过率。i(

1.2　分裂窗算法各中间参数的计算1.2.1　亮度温度的计算

MODIS图像是用DN值表示的,因此,要计算星上亮

式中,band2和band19分别为第2和19波段的DN值;

scale2、offset2、scale19和offset19是定标常量,可从MODIS

科学数据集的属性数据中查出。

大气透过率除受大气水汽含量影响外,还受遥感器视角和大气剖面温度的影响,因此大气透过率的估计还需要进行视角和温度校正

(7)(8)

[5～6]

温,必须先将DN值转换成相应的辐射强度值,然后再用

Planck函数求解星上亮温。

MODIS第31～32波段的辐射强度值计算公式如下:

rad31=scale32(band31-offset31)rad32=scale32(band32-offset32)

。MODIS第31和32波段

(14)

大气透过率可用下式来计算:

τθ)=τττθ)i(i(10)+δi(T)-δi(

θ)是第i(i=31,32)波段的大气透过率;τ式中,τi(i(10)为传感器视角为10°时的星下大气透过率,可根据表1给τ出的方程由大气水汽含量计算得出;δi(T)是大气透过τθ)率的温度校正函数,由表2给出的方程计算得出;δi(是传感器视角校正函数,由公式(18)和公式(19)给出。

式中,rad31和rad32分别为MODIS第31～32波段的热辐

-2-1-1

射强度(Wmsr(m);band31、band32分别为MODIS

第31、32波段的DN值;scale31和offset31为MODIS第31波段的辐射定标常量,可从MODIS数据集的属性数据中

表1　MODIS第31和32波段的星下大气透过率估计方程

Tab.1　Theestimationequationofatmospherictransmittanceusingband31and32ofMODIS波段

MODIS31MODIS32

水汽含量0.4～2.0(g/cm2)

τ31(10)=0.99513-0.08082wτ32(10)=0.99376-0.11369w

水汽含量2.0～4.0(g/cm2)

τ31(10)=1.08692-0.12759wτ32(10)=1.07900-0.15925w

水汽含量4.0～6.0(g/cm2)

τ31(10)=1.07268-0.12571wτ32(10)=0.93821-0.12613w

表2　大气透过率的温度校正函数

Tab.2　Thetemperaturerevisedfunctionofatmospherictransmittance

波段

MODIS31MODIS32

T>318K

278

δτ31(T)=0.08δτ32(T)=0.095

δτ31(T)=-0.05+0.00325(T31-278)δτ32(T)=-0.065+0.004(T32-278)

δτ31(T)=-0.05δτ32(T)=-0.065

　　　　注:T是第31和32波段的亮度温度。

大气透过率的传感器视角校正函数如下:

-52

δτθ)=-0.00322+(3.0967×10)θi(

HDF_SD_ATTRINFO,sds_id,8,data=scale;　读取

(15)

scale属性数据。

HDF_SD_ATTRINFO,sds_id,9,data=offset;　读取offset属性数据。

HDF_SD_GETDATA,sds_id,data;　读取科学数据集

式中,θ是传感器视角,可由sensorzenth数据集提供。

1.2.3　地表比辐射率的计算

地表比辐射率的计算要分陆地像元和水体像元两种情况。对于水体像元可直接取ε.99683,ε31=032=

0.99254。对于陆地像元地表比辐射率可由植被覆盖率计

各波段的数据。

band1=data[3,3,0];　把第1波段数据存入band1中。

HDF_SD_ENDACCESS,sds_id;　关闭SDdatasetHDF_SD_END,sd_id;　关闭文件2.2　大气透过率的算法实现

算得出。因此,要计算陆地像元地表比辐射率,必须先求出植被覆盖率。植被覆盖率的求解公式如下

[7]

:

(16)

pv=(NDVI-NDVIs)/(NDVIv-NDVIs)

式中,NDVI是植被指数,NDVIv和NDVIs分别是茂密植被覆盖和完全裸土像元的NDVI值,通常取NDVIv=0.9,

NDVIs=0.15。对于MODIS图像而言,NDVI可用第1和2

,根据

1.2波段来计算:

NDVI=(ref2-ref1)/(ref2-ref1)

(17)

)

。

式中,ref1和ref2分别是MODIS图像第1和率。

6]

:

(18)(19)

εε31=vvv+(1-Pv)Rsε31s+dεε32=PvRvε32v+(1-Pv)Rsε32s+d

式中,εS第31,32波段的地表比辐射率;31和ε32是MODIε31v和ε31s分别是植被和裸土在第31波段的地表比辐射率,分别取ε.98672,ε.96767;ε31v=031s=032v和ε32s分别是植被和裸土在第32波段的地表比辐射率,分别取ε32v=

0.98990,ε.97790;Pv是像元的植被覆盖率;Rv和32s=0

Rs分别是植被和裸土的辐射比率,可取Rv=0192762+

图1　大气透过率计算流程图

Fig.1　Thecalculationworkflowchartfor

atmospherictransmittance

ε是热辐射相互0.07033Pv,Rs=0.99782+0.08362Pv;d

作用校正,由植被和裸土之间的热辐射相互作用产生,可通过如下经验公式来估计

[7]

计算像元的星下大气透过率时,首先要判断该像元的大气水汽含量处于哪个范围,然后再根据相应的公式计算。这可以通过循环语句和if语句来实现,但是IDL提供的where()函数能更方便完成此任务,而且程序的运算速度更快。where()函数返回数组或数组表达式中非零元素的下标。在IDL中任何数组都可看作是1维数组,因此当where()函数的自变量是2维或更高维数组时,返回的带有下标的数组总是1维的,即将输入数组看作1维数组。我们可以利用where函数查找water_vapor数组中水汽含量符合某条件的元素,返回它们的下标,然后再利用返回的下标数组,对星下大气透过率数组中相应元素进行赋值。

2.3　地表比辐射率的算法实现

:

(20)

ε=0.003796min(Pv,(1-Pv))d

式中,min(Pv,(1-Pv))表示取Pv和1-Pv的最小值。

2　算法实现与优化

2.1　MODIS定标数据与DN值的读取

用分裂窗算法反演地表温度所需数据为1、2、19、31、

32波段的定标数据和DN值。IDL提供的HDF_SD_AT2TRINFO过程可用于读取HDF科学数据集的属性数据,HDF_SD_GETDATA过程可用于读取HDF科学数据集的

数据。下面以第1波段为例,介绍读取定标数据和DN值的代码(分号后为注释文字)。

sd_id=HDF_SD_START(filename,/read);　打开文

　　地表比辐射率是分裂窗算法的又一重要参数,它的计算需要区分水体像元和陆地像元。我们采用的算法是:首先用MODIS第1、2波段计算NDVI,然后用where()函数来查找NDVI

NDVI>0的像元(陆地像元),计算这些像元的植被覆盖

件,并初始化SD界面

sds_id=HDF_SD_SELECT(sd_id,pos);　pos为科学

数据集索引值(如表3所示)。

HDF_SD_GETINFO,sds_id,dims=dims;　获取图像

的维度信息。

率pv,进而计算其地表比辐射率。算法流程如图2

所示。

存函数IDL_MemFree()实现,也可以通过另一个更为简便的函数temporary(),这个函数在进行两个或者多个数组计算时,只是对数组自身进行一个拷贝,无需产生真实的需要分配内存的中间变量,所以能节省大量的空间和运算时间。

3结　论

本文基于覃志豪提出的适用于MODIS数据的地表温

度反演分裂窗算法,在IDL6.0环境下,编程实现了MO2

图2　地表比辐射率计算流程图

Fig.2　Thecalculationworkflowchartfor

landsurfaceradiance

DIS数据地表温度反演的全自动化运算。结果表明,本文

提出的MODIS数据地表温度反演的实现方法,操作简便、运算速度快,可用于批量、快速、MODIS地表温度。

2.4　地表温度的算法实现与结果保存

计算出大气透过率和地表比辐射率这两个参数后,地表温度就可用分裂窗算法(公式(1)-(8))IDLL中的[　OGD,KarnieliA,eta1.Derivationofsplitwindow

algorithmanditssensitivityanalysisforretrievinglandsurfacetem2peraterefromNOAA-AVHRRdata[J].JournalofGeophysicalResearch,2001,106(D19):22655-22670.

[2]　覃志豪,高懋芳,秦晓敏,等.农业旱灾监测中的地表温度遥

读写函数,:

SD_id=(me,/CREATE);　

创建一个HDFSDS_id=HDF_SD_CREATE(sd_id,′LST′,[dims[0],dims[1]],/FLOAT);　创建并定义HDF文件的科学数据

感反演方法———以MODIS数据为例[J].自然灾害学报,

2005,14(3).

[3]　毛克彪,覃志豪.MODIS影像反演环渤海地区的大气水汽含

集。

HDF_SD_SETINFO,sds_id,FILL=0.0,LABEL=′LST′,RANGE=[dims[0],dims[1]];　向创建的科学数

量[J].遥感信息,2004,76(4).

[4]　KaufmanYJ,andGaoBC.Remotesensingofwatervaporinthe

nearIRfromEOS/MODIS[J].IEEETransactionsonGeosciencesandRemoteSensing,30(5):871-884.

[5]　KerrYH,LagouardeJP,andImbernonJ.Accuratelandsurface

temperatureretrievalfromAVHRRdatawithuseofanimprovedsplitwindowalgorithm[J].1992,41:197-209.

[6]　覃志豪,李文娟,吴波,等.LandsatTM6波段范围内地表比辐

RemoteSensingofEnvironment,

据集写入相关信息。

HDF_SD_ADDDATA,sds_id,LST;　向科学数据集添

加LST数据

HDF_SD_ENDACCESS,sds_id;　关闭科学数据集HDF_SD_END,sd_id;　关闭文件2.5　算法的优化处理

由于MODIS数据量比较大,分裂窗算法又比较复杂,所需中间变量多,而IDL采用的是面向数组的程序运行方式,数据量多大就要开辟多大的内存来存储这些数据,因此,程序在运行过程中,内存开销将会很大。为了节省计算机的系统资源,提高程序的运行效率和速度,必须对程序进行优化处理。我们采取的具体措施如下:

1)仅读取1、2、19、31、32波段数据,并且是用到哪个

射率的估计[J].国土资源遥感,2004,(3):27-32.

[7]　SobrinoJA,Jimenez-MuiiozJC,andaoliniL.Landsurface

temperatureretrievalfromLANDSATTM5[J].RemoteSensingofEnvironment,2004,90:434-440.

[8]　WanZ,DozierJ.Ageneralizedsplit-windowsalgorithmforre2

trievinglandsurfacetemperaturefromspace[J].

IEEETransac2

tionsonGeosciencesandRemoteSensing,1996,34:892-905.[9]　sobrinoJA,CollC,andCasellesV.Atmosphericcorrectionfor

landsurfacetemperatureusingNOAA-11AVHRRchannels4and5[J].RemoteSensingofEnvironment,1991,38:19-34.

波段的数据再进行读取,而不是程序开始时就把所有波段数据都读进来。

2)在程序运算中尽量收回中间变量所占的内存空

间,这可以通过分配内存函数IDL_MemAlloc()和释放内

[责任编辑:栾丽杰]